Разрешаем совместную работу на Python и R:
%load_ext rpy2.ipython
Устанавливаем необходимые библиотеки Python:
import pip
import warnings
warnings.simplefilter('ignore')
warnings.filterwarnings('ignore') # предупреждающие сообщения не выводятся
pip.main(['install', 'pandas'])
import pandas as pd
pip.main(['install', 'matplotlib'])
import matplotlib
import matplotlib.pyplot as plt
#графики отображаются сразу после кода
%matplotlib inline
pip.main(['install', 'numpy'])
import numpy as np
pip.main(['install', 'seaborn'])
import seaborn as sns
sns.set(style='darkgrid')
pip.main(['install', 'math'])
import math
pip.main(['install', 'ast'])
import ast # literal_eval
pip.main(['install', 'outlier_utils'])
from outliers import smirnov_grubbs as grubbs # критерий Граббса
pip.main(['install', 'statsmodels'])
import statsmodels.api as sm # qq-график
from statsmodels.stats.diagnostic import lilliefors # критерий Колмогорова-Смирнова в мод. Лиллиефорса
from statsmodels.stats.contingency_tables import mcnemar # тест МакНемара
from statsmodels.stats.outliers_influence import variance_inflation_factor # фактор инфляции дисперсии
from statsmodels.tools.tools import add_constant
from statsmodels.formula.api import ols # многофакторный дисперсионный анализ
pip.main(['install', 'scipy'])
import scipy.stats as stats
from scipy.stats import kstest # критерий Колмогорова-Смирнова
from scipy.stats import shapiro # критерий Шапиро_Уилка
from scipy.stats import anderson # критерий Андерсона-Дарлинга
from scipy.stats import anderson_ksamp
from scipy.stats import cramervonmises # критерий Крамера фон Мизеса
from scipy.stats import cramervonmises_2samp
from scipy.stats import mannwhitneyu # критерий Уилкоксона-Манна-Уитни
from scipy.stats import levene # критерий Левене
from scipy.stats import bartlett # критерий Баретта
from scipy.stats import fligner # критерий Флигнера-Килина
from scipy.stats import pearsonr # коэффициент корреляции Пирсона
from scipy.stats import spearmanr # коэффициент корреляции Спирмена
from scipy.stats import kendalltau # коэффициент корреляции Кендалла
from scipy.stats import chi2_contingency # критерий хи-квадрат
from scipy.stats import fisher_exact # точный тест Фишера
from scipy.stats import f_oneway # однофакторный дисперсионный анализ
pip.main(['install', 'sfrancia'])
from sfrancia import shapiroFrancia # критерий Шапиро-Франсия
pip.main(['install', 'statistics'])
import statistics # mean и т.д.
pip.main(['install', 'pingouin'])
from pingouin import ttest # критерий Стьюдента
pip.main(['install', 'cmh'])
from cmh import CMH # критерий Кохрана-Мантеля-Хензеля
/usr/local/lib/python3.10/dist-packages/_distutils_hack/__init__.py:33: UserWarning: Setuptools is replacing distutils.
warnings.warn("Setuptools is replacing distutils.")
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: pandas in /usr/local/lib/python3.10/dist-packages (1.5.3)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas) (2023.3.post1)
Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas) (1.23.5)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.8.1->pandas) (1.16.0)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: matplotlib in /usr/local/lib/python3.10/dist-packages (3.7.1)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (4.45.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.4.5)
Requirement already satisfied: numpy>=1.20 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.23.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (23.2)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (3.1.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (2.8.2)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (1.23.5)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: seaborn in /usr/local/lib/python3.10/dist-packages (0.12.2)
Requirement already satisfied: numpy!=1.24.0,>=1.17 in /usr/local/lib/python3.10/dist-packages (from seaborn) (1.23.5)
Requirement already satisfied: pandas>=0.25 in /usr/local/lib/python3.10/dist-packages (from seaborn) (1.5.3)
Requirement already satisfied: matplotlib!=3.6.1,>=3.1 in /usr/local/lib/python3.10/dist-packages (from seaborn) (3.7.1)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (4.45.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (23.2)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (3.1.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=0.25->seaborn) (2023.3.post1)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib!=3.6.1,>=3.1->seaborn) (1.16.0)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
ERROR: Could not find a version that satisfies the requirement math (from versions: none)
ERROR: No matching distribution found for math
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Collecting ast
Using cached AST-0.0.2.tar.gz (19 kB)
Preparing metadata (setup.py): started
error: subprocess-exited-with-error × python setup.py egg_info did not run successfully. │ exit code: 1 ╰─> See above for output. note: This error originates from a subprocess, and is likely not a problem with pip.
Preparing metadata (setup.py): finished with status 'error'
error: metadata-generation-failed × Encountered error while generating package metadata. ╰─> See above for output. note: This is an issue with the package mentioned above, not pip. hint: See above for details.
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: outlier_utils in /usr/local/lib/python3.10/dist-packages (0.0.5)
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from outlier_utils) (1.23.5)
Requirement already satisfied: scipy in /usr/local/lib/python3.10/dist-packages (from outlier_utils) (1.11.4)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: statsmodels in /usr/local/lib/python3.10/dist-packages (0.14.0)
Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (1.23.5)
Requirement already satisfied: scipy!=1.9.2,>=1.4 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (1.11.4)
Requirement already satisfied: pandas>=1.0 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (1.5.3)
Requirement already satisfied: patsy>=0.5.2 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (0.5.3)
Requirement already satisfied: packaging>=21.3 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (23.2)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.0->statsmodels) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.0->statsmodels) (2023.3.post1)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from patsy>=0.5.2->statsmodels) (1.16.0)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: scipy in /usr/local/lib/python3.10/dist-packages (1.11.4)
Requirement already satisfied: numpy<1.28.0,>=1.21.6 in /usr/local/lib/python3.10/dist-packages (from scipy) (1.23.5)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: sfrancia in /usr/local/lib/python3.10/dist-packages (1.0.8)
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from sfrancia) (1.23.5)
Requirement already satisfied: scipy in /usr/local/lib/python3.10/dist-packages (from sfrancia) (1.11.4)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: statistics in /usr/local/lib/python3.10/dist-packages (1.0.3.5)
Requirement already satisfied: docutils>=0.3 in /usr/local/lib/python3.10/dist-packages (from statistics) (0.18.1)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: pingouin in /usr/local/lib/python3.10/dist-packages (0.5.3)
Requirement already satisfied: numpy>=1.19 in /usr/local/lib/python3.10/dist-packages (from pingouin) (1.23.5)
Requirement already satisfied: scipy>=1.7 in /usr/local/lib/python3.10/dist-packages (from pingouin) (1.11.4)
Requirement already satisfied: pandas>=1.0 in /usr/local/lib/python3.10/dist-packages (from pingouin) (1.5.3)
Requirement already satisfied: matplotlib>=3.0.2 in /usr/local/lib/python3.10/dist-packages (from pingouin) (3.7.1)
Requirement already satisfied: seaborn>=0.11 in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.12.2)
Requirement already satisfied: statsmodels>=0.13 in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.14.0)
Requirement already satisfied: scikit-learn in /usr/local/lib/python3.10/dist-packages (from pingouin) (1.2.2)
Requirement already satisfied: pandas-flavor>=0.2.0 in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.6.0)
Requirement already satisfied: outdated in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.2.2)
Requirement already satisfied: tabulate in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.9.0)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (4.45.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (23.2)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (3.1.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.0->pingouin) (2023.3.post1)
Requirement already satisfied: xarray in /usr/local/lib/python3.10/dist-packages (from pandas-flavor>=0.2.0->pingouin) (2023.7.0)
Requirement already satisfied: patsy>=0.5.2 in /usr/local/lib/python3.10/dist-packages (from statsmodels>=0.13->pingouin) (0.5.3)
Requirement already satisfied: setuptools>=44 in /usr/local/lib/python3.10/dist-packages (from outdated->pingouin) (67.7.2)
Requirement already satisfied: littleutils in /usr/local/lib/python3.10/dist-packages (from outdated->pingouin) (0.2.2)
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (from outdated->pingouin) (2.31.0)
Requirement already satisfied: joblib>=1.1.1 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->pingouin) (1.3.2)
Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->pingouin) (3.2.0)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from patsy>=0.5.2->statsmodels>=0.13->pingouin) (1.16.0)
Requirement already satisfied: charset-normalizer<4,>=2 in /usr/local/lib/python3.10/dist-packages (from requests->outdated->pingouin) (3.3.2)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests->outdated->pingouin) (3.6)
Requirement already satisfied: urllib3<3,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests->outdated->pingouin) (2.0.7)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests->outdated->pingouin) (2023.11.17)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip. Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue. To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: cmh in /usr/local/lib/python3.10/dist-packages (1.0.1)
Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.10/dist-packages (from cmh) (1.23.5)
Requirement already satisfied: pandas>=1 in /usr/local/lib/python3.10/dist-packages (from cmh) (1.5.3)
Requirement already satisfied: scikit-learn in /usr/local/lib/python3.10/dist-packages (from cmh) (1.2.2)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1->cmh) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1->cmh) (2023.3.post1)
Requirement already satisfied: scipy>=1.3.2 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->cmh) (1.11.4)
Requirement already satisfied: joblib>=1.1.1 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->cmh) (1.3.2)
Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->cmh) (3.2.0)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.8.1->pandas>=1->cmh) (1.16.0)
Скачиваем функции, которые понадобятся для работы на R:
%%R
oldw <- getOption("warn") # предупреждающие сообщения не выводятся
options(warn = -1)
install.packages("outliers") # поиск выбросов (критерии Граббса и Диксона)
library(outliers)
install.packages("nortest") # критерии Андерсона-Дарлинга, Лиллиефорса, Шапиро-Франсия; Пирсона
library(nortest)
install.packages("goftest") # критерий Крамера фон Мизеса
library(goftest)
install.packages("stats") # rnorm, qq-график
library(stats)
install.packages("boot") # метод огибающих
library(boot)
install.packages("corrplot") # корреляционная матрица
library(corrplot)
install.packages("car") # vif
library(car)
install.packages("data.table") # %like%
library(data.table)
install.packages("dplyr") # %>%
library(dplyr)
install.packages("relevance") # dropNA
library(relevance)
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/outliers_0.15.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 15967 bytes (15 KB)
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 15 KB
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/nortest_1.0-4.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 6179 bytes
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 6179 bytes
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/goftest_1.2-3.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 32319 bytes (31 KB)
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 31 KB
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]:
Attaching package: ‘goftest’
WARNING: R[write to console]: The following objects are masked from ‘package:nortest’:
ad.test, cvm.test
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/boot_1.3-28.1.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 236854 bytes (231 KB)
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 231 KB
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/corrplot_0.92.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 3765850 bytes (3.6 MB)
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 3.6 MB
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: corrplot 0.92 loaded
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/car_3.1-2.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 579829 bytes (566 KB)
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 566 KB
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: Loading required package: carData
WARNING: R[write to console]:
Attaching package: ‘car’
WARNING: R[write to console]: The following object is masked from ‘package:boot’:
logit
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/data.table_1.14.10.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 5311651 bytes (5.1 MB)
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 5.1 MB
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: data.table 1.14.10 using 1 threads (see ?getDTthreads). Latest news: r-datatable.com
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/dplyr_1.1.4.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 1207521 bytes (1.2 MB)
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 1.2 MB
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]:
Attaching package: ‘dplyr’
WARNING: R[write to console]: The following objects are masked from ‘package:data.table’:
between, first, last
WARNING: R[write to console]: The following object is masked from ‘package:car’:
recode
WARNING: R[write to console]: The following objects are masked from ‘package:stats’:
filter, lag
WARNING: R[write to console]: The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/relevance_2.0.tar.gz'
WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]: length 889295 bytes (868 KB)
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]:
WARNING: R[write to console]: downloaded 868 KB
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]: The downloaded source packages are in
‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]:
WARNING: R[write to console]:
WARNING: R[write to console]:
Attaching package: ‘relevance’
WARNING: R[write to console]: The following object is masked from ‘package:dplyr’:
last
WARNING: R[write to console]: The following object is masked from ‘package:data.table’:
last
В работе (в различных заданиях) использованы 2 датасета.
Первый набор данных основан на датасете best-selling-books. В нём содержится информация о самых продаваемых книгах (бестселлерах[1](#01)):
Во-первых, для удобства обработки произведения, написанные в соавторстве, разделены на строки, соответствующие отдельным авторам.
Во-вторых, данные были дополнены следующей информацией, взятой из открытых источников:
В-третьих, из-за большого количества пропусков и запутанных обозначений был полностью изменён столбец genres. Для единообразия обозначений (в том числе со вторым датасетом) использовались данные сайта Goodreads.
1: Под ними понимаются те книги, которые имеют наибольшее ожидаемое количество проданных экземпляров, а не количество экземпляров, напечатанных или находящихся в собственности в настоящее время; также в список не входят комиксы и учебная литература.
Получившийся набор данных выглядит следующим образом:
Python:
df1 = pd.read_csv('best-selling-books-upgrade.csv')
# заполняем пропуски в genres
df1.genres = df1.genres.fillna('[]')
# преобразуем данные столбца genres из строк в списки
df1.genres = df1.genres.apply(ast.literal_eval)
print('Количество строк:', df1.shape[0])
df1.head(5)
Количество строк: 175
| title | author | gender | birth_year | age | language | first_published | sales | genres | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | And Then There Were None | Agatha Christie | F | 1890 | 49 | English | 1939 | 100.0 | [Classics, Fiction, Crime, Thriller, Mystery T... |
| 1 | The Plague (La Peste) | Albert Camus | M | 1913 | 34 | French | 1947 | 12.0 | [Fiction, Classics, Philosophy, France, Litera... |
| 2 | The Stranger (L'Étranger) | Albert Camus | M | 1913 | 29 | French | 1942 | 10.0 | [Classics, Fiction, Philosophy, France, Litera... |
| 3 | The Joy of Sex | Alex Comfort | M | 1920 | 52 | English | 1972 | 10.0 | [Nonfiction, Sexuality, Relationships, Self He... |
| 4 | The Young Guard (Молодая гвардия) | Alexander Alexandrovich Fadeyev | M | 1901 | 44 | Russian | 1945 | 26.0 | [Russia, Fiction, Russian Literature, Historic... |
R:
%%R
df1 <- read.csv(file = "best-selling-books-upgrade.csv", header = TRUE)
# заполняем пропуски в genres
df1$genres[is.na(df1$genres)] <- "[]"
# преобразуем данные столбца genres из строк в списки
df1$genres <- gsub("\\[|\\]|\\'|\\ ", "", df1$genres)
df1$genres <- strsplit(df1$genres, "\\,")
cat("Количество строк:", nrow(df1), "\n")
head(df1, n=5)
Количество строк: 175
title author gender
1 And Then There Were None Agatha Christie F
2 The Plague (La Peste) Albert Camus M
3 The Stranger (L'Étranger) Albert Camus M
4 The Joy of Sex Alex Comfort M
5 The Young Guard (Молодая гвардия) Alexander Alexandrovich Fadeyev M
birth_year age language first_published sales
1 1890 49 English 1939 100
2 1913 34 French 1947 12
3 1913 29 French 1942 10
4 1920 52 English 1972 10
5 1901 44 Russian 1945 26
genres
1 Classics, Fiction, Crime, Thriller, MysteryThriller, Audiobook, Mystery
2 Fiction, Classics, Philosophy, France, Literature, FrenchLiterature, Novels
3 Classics, Fiction, Philosophy, France, Literature, Novels, FrenchLiterature
4 Nonfiction, Sexuality, Relationships, SelfHelp, Reference, Erotica, Psychology
5 Russia, Fiction, RussianLiterature, HistoricalFiction, War, Classics, WorldWarII, Novels, Unfinished, Childrens
Второй использованный датасет - Best Books (10k) Multi-Genre Data. В нём содержится информация о 10,000 самых рекоммендованных книгах согласно рейтингу сайта Goodreads. В работе используются следующие данные:
Так выглядит этот набор данных:
Python:
df2 = pd.read_csv('goodreads-data-upgrade.csv')
# заполняем пропуски в genres
df2.genres = df2.genres.fillna('[]')
# преобразуем данные столбца genres из строк в списки
df2.genres = df2.genres.apply(ast.literal_eval)
# преобразуем данные столбца number_of_ratings из строк в числа
df2['number_of_ratings'] = df2['number_of_ratings'].apply(lambda x: int(x.replace(',', '')))
print('Количество строк:', df2.shape[0])
df2.head(5)
Количество строк: 10000
| title | author | average_rating | number_of_ratings | genres | |
|---|---|---|---|---|---|
| 0 | To Kill a Mockingbird | Harper Lee | 4.27 | 5691311 | [Classics, Fiction, Historical Fiction, School... |
| 1 | Harry Potter and the Philosopher’s Stone (Harr... | J.K. Rowling | 4.47 | 9278135 | [Fantasy, Fiction, Young Adult, Magic, Childre... |
| 2 | Pride and Prejudice | Jane Austen | 4.28 | 3944155 | [Classics, Fiction, Romance, Historical Fictio... |
| 3 | The Diary of a Young Girl | Anne Frank | 4.18 | 3488438 | [Classics, Nonfiction, History, Biography, Mem... |
| 4 | Animal Farm | George Orwell | 3.98 | 3575172 | [Classics, Fiction, Dystopia, Fantasy, Politic... |
R:
%%R
df2 <- read.csv(file = "goodreads-data-upgrade.csv", header = TRUE)
# заполняем пропуски в genres
df2$genres[is.na(df2$genres)] <- "[]"
# преобразуем данные столбца genres из строк в списки
df2$genres <- gsub("\\[|\\]|\\'|\\ ", "", df2$genres)
df2$genres <- strsplit(df2$genres, "\\,")
# преобразуем данные столбца number_of_ratings из строк в числа
df2$number_of_ratings <- as.integer(gsub(",", "", df2$number_of_ratings))
cat("Количество строк:", nrow(df2), "\n")
head(df2, n=5)
Количество строк: 10000
title author
1 To Kill a Mockingbird Harper Lee
2 Harry Potter and the Philosopher’s Stone (Harry Potter, #1) J.K. Rowling
3 Pride and Prejudice Jane Austen
4 The Diary of a Young Girl Anne Frank
5 Animal Farm George Orwell
average_rating number_of_ratings
1 4.27 5691311
2 4.47 9278135
3 4.28 3944155
4 4.18 3488438
5 3.98 3575172
genres
1 Classics, Fiction, HistoricalFiction, School, Literature, YoungAdult, Historical
2 Fantasy, Fiction, YoungAdult, Magic, Childrens, MiddleGrade, Classics
3 Classics, Fiction, Romance, HistoricalFiction, Literature, Historical, Audiobook
4 Classics, Nonfiction, History, Biography, Memoir, Historical, Holocaust
5 Classics, Fiction, Dystopia, Fantasy, Politics, School, Literature
Ядерные оценки позволяют строить гладкие оценки плотностей распределения. Часто их изображают вместе с соответствующими гистограммами. Рассмотрим, например, распределение годов первой публикации бестселлеров.
Для определения количества интервалов разбиения hist() в R по умолчанию использует метод Стёрджеса , а histplot() в Python - максимум оценок Стёрджеса и Фридмана-Диакониса. Метод Стёрджеса в данном случае даёт очевидно неудовлетворительный результат, поэтому в обоих случаях использован метод Фридмана-Диакониса.
Тем не менее, отличия в графиках присутствуют. Видимо, метод по-разному реализован на R и Python: границы отрезков разбиения и их количество не совпадают.
Полоса пропускания ядра установлена по умолчанию.
Python:
edges = np.histogram_bin_edges(df1['first_published'], bins='fd')
print('Количество отрезков разбиения:', edges.shape[0] - 1)
edges
Количество отрезков разбиения: 42
array([1304., 1321., 1338., 1355., 1372., 1389., 1406., 1423., 1440.,
1457., 1474., 1491., 1508., 1525., 1542., 1559., 1576., 1593.,
1610., 1627., 1644., 1661., 1678., 1695., 1712., 1729., 1746.,
1763., 1780., 1797., 1814., 1831., 1848., 1865., 1882., 1899.,
1916., 1933., 1950., 1967., 1984., 2001., 2018.])
plt.figure(figsize=(10, 6))
sns.histplot(data=df1['first_published'], stat='density', bins=edges, kde=True)
plt.xticks(fontsize=9)
plt.xlabel('Год', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('Распределение годов первой публикации', fontsize=10, fontweight='bold')
plt.show()
R:
%%R
edges = hist(df1$first_published, freq = FALSE, breaks = "FD", plot = FALSE)
cat("Количество отрезков разбиения:", length(edges$breaks) - 1, "\n")
edges$breaks
Количество отрезков разбиения: 36 [1] 1300 1320 1340 1360 1380 1400 1420 1440 1460 1480 1500 1520 1540 1560 1580 [16] 1600 1620 1640 1660 1680 1700 1720 1740 1760 1780 1800 1820 1840 1860 1880 [31] 1900 1920 1940 1960 1980 2000 2020
%%R -w 1000 -h 600
hist(df1$first_published, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Год", ylab = "",
main = "Распределение годов первой публикации")
lines(density(df1$first_published, adjust=), col = "royalblue", lwd = 2)
Видно, что пик распределения приходится примерно на 80е - 90е годы. Действительно это время, когда, с одной стороны, книги активно пишутся и издаются и уже давно доступны большому количеству людей, но, с другой стороны, ещё не распространён интернет и т.д.
Если строить ядерные оценки плотностей для разных групп, то можно сравнивать распределения в них. Рассмортим, например, как отличаются распределения для авторов разного пола.
Признак gender выбран намеренно, так как выделенные с помощью него подгруппы не пересекаются и любая книга обязательно попадает в одну из подгрупп. В противном случее использовать нормирование (площадь под гистограммой равна 1) было бы неразумно (например, из-за того, что большинство книг имеют не один жанр, некоторые из них при построении отдельных графиков могут учитываться несколько раз и т.д.).
Интервалы разбиения у графиков такие же, как и у графиков выше, полоса пропускания по умолчанию.
Python:
plt.figure(figsize=(10, 6))
sns.histplot(data=df1[df1.gender == 'F'].first_published, stat='density', kde=True,
bins=edges,
color='red', alpha=0.3,
label='Female')
sns.histplot(data=df1[df1.gender == 'M'].first_published, stat='density', kde=True,
bins=edges,
color='blue', alpha=0.3,
label='Male')
plt.xticks(fontsize=9)
plt.xlabel('Год', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('Распределение годов первой публикации', fontsize=10, fontweight='bold')
plt.legend(title='Пол', title_fontsize=10, fontsize=10)
plt.show()
R:
%%R
# функция, которая создаёт прозрачные версии цветов
t_col <- function(color, percent = 50, name = NULL) {
# color = color name
# percent = % transparency
# name = an optional name for the color
# Get RGB values for named color
rgb.val <- col2rgb(color)
# Make new color using input color as base and alpha set by transparency
t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3],
max = 255,
alpha = (100 - percent) * 255 / 100,
names = name)
# Save the color
invisible(t.col)
}
%%R -w 1000 -h 600
breaks = hist(df1$first_published, breaks = 145, plot = FALSE)$breaks
my_red <- t_col("red", perc = 60)
hist(df1[df1$gender %like% "F", "first_published"], freq = FALSE,
breaks = edges$breaks,
xlim = c(1300, 2020),
col = my_red,
xlab = "Год", ylab = "",
main = "Распределение годов первой публикации")
my_blue <- t_col("royalblue", perc = 60)
hist(df1[df1$gender %like% "M", "first_published"], freq = FALSE,
breaks = edges$breaks,
col = my_blue,
add = TRUE)
lines(density(df1[df1$gender %like% "F", "first_published"], adjust=), col = "red", lwd = 3)
lines(density(df1[df1$gender %like% "M", "first_published"], adjust=), col = "royalblue", lwd = 3)
legend("topleft", inset=c(0.03, 0), legend=c("Female", "Male"), fill = c(my_red, my_blue), title="Пол")
На графиках видно, что со временем растёт количество писательниц.
Для визуализации связи между двумя переменными, одна из которых является количественной, а другая качественной, можно использовать график условной плотности. Рассмотрим, например, какова связь между языком и количеством проданных копий.
Так как <<оценки условных плотностей более надёжны для областей с высокой плотностью $x$>>, а в данных наблюдается <<перекос>> (см. stripchart), то отображать на графике будем только часть данных - исключим выбросы справа (см. boxplot). В противном случае на части, где данных мало, график начинает выглядеть неадекватно.
Python:
lang_colors = reversed(["#B81840", "#C43142", "#D04544", "#DB5745", "#E56946",
"#EE7D47", "#F2994E", "#F5B355", "#F3A469", "#EC867E",
"#DD708E", "#BD6C99", "#9A68A4", "#76669B", "#4D6291", "#26456E"])
plt.figure(figsize=(10, 6))
ax = sns.kdeplot(data=df1[df1.sales < 66], x='sales', hue='language',
multiple='fill',
palette=lang_colors)
plt.xticks(fontsize=9)
plt.xlabel('Количество проданных экземпляров', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('Распределение количества проданных экземпляров по языкам', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.9), borderaxespad=0,
title='Язык', title_fontsize=10, fontsize=10)
ax.set_xlim(min(df1.sales), 66)
plt.show()
R:
%%R
lang_colors <- c("#B81840", "#C43142", "#D04544", "#DB5745", "#E56946",
"#EE7D47", "#F2994E", "#F5B355", "#F3A469", "#EC867E",
"#DD708E", "#BD6C99", "#9A68A4", "#76669B", "#4D6291", "#26456E")
%%R -w 1000 -h 600
# добавляем дополнительное место слева от графика для легенды
par(mar=c(5, 4, 4, 16), xpd=TRUE)
languages <- factor(df1[df1$sales < 66, ]$language)
lang_levels <- unique(languages)
sales <- df1[df1$sales < 66, ]$sales
cdplot(languages ~ sales, col = lang_colors, ylevels = lang_levels,
xlab = "Количество проденных экземпляров", ylab="", yaxlab = "",
main = "Распределение количества проданных экземпляров по языкам")
# добавляем метки осей слева
#axis(2)
# добавляем легенду вне графика
legend("topright", inset=c(-0.2, 0.25), legend=lang_levels, fill = rev(lang_colors), title="Язык")
На графиках хорошо видно, что книги на английском языке чаще становятся бестселлерами.
Одна из альтернатив столбчатым диаграммам - точечные диаграммы Кливленда. С их помощью рассмотрим, как связаны продажи книги с языком, на котором она написана. Для каждого языка найдём книгу с максимальным количеством проданных экземпляров и нанесём соответствующее значение на график (причём отсортируем языки в порядке возрастания количества произведений, написанных на том же языке).
Python:
language_freq = df1['language'].value_counts().reset_index()
language_freq.columns = ['language', 'freq']
language_max_sales = df1.groupby('language')['sales'].max().reset_index()
language_max_sales.columns = ['language', 'max_sales']
tmp1_df = pd.merge(language_freq, language_max_sales)
tmp1_df = tmp1_df.sort_values(['freq', 'max_sales'], ascending = (True, True))
plt.figure(figsize=(10, 6))
ax = sns.stripplot(x=tmp1_df.max_sales, y=tmp1_df.language, hue=tmp1_df.freq, palette='deep', size=7)
plt.xticks(fontsize=9)
plt.xlabel('Максимальное количество проданных экземпляров, млн', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Язык', fontsize=10)
plt.title('Связь максимального количества проданных экземпляров с количеством произведений на языке', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0, title='Количество произведений', fontsize=10, title_fontsize=10)
plt.show()
R:
%%R
language_freq <- as.data.frame(table(df1$language))
names(language_freq) <- c("language", "freq")
language_max_sales <- df1[order(df1$language, - df1$sales), ]
language_max_sales <- language_max_sales[!duplicated(language_max_sales$language), ]
language_max_sales <- language_max_sales[ , c("language", "sales")]
tmp1_df <- merge(language_max_sales, language_freq, by="language")
tmp1_df <- tmp1_df[order(-tmp1_df$freq, -tmp1_df$sales), ]
names(tmp1_df) <- c("language", "max_sales", "freq")
tmp1_df[tmp1_df$freq == 1, "color"] <- "royalblue"
tmp1_df[tmp1_df$freq == 2, "color"] <- "yellow"
tmp1_df[tmp1_df$freq == 3, "color"] <- "green"
tmp1_df[tmp1_df$freq == 4, "color"] <- "red"
tmp1_df[tmp1_df$freq == 5, "color"] <- "purple"
tmp1_df[tmp1_df$freq == 6, "color"] <- "brown"
tmp1_df[tmp1_df$freq == 132, "color"] <- "pink"
%%R -w 1000 -h 600
dotchart(tmp1_df$max_sales,
labels = tmp1_df$language, groups = as.factor(tmp1_df$freq),
pch = 16, cex = 1.2,
color = tmp1_df$color, gcolor = "black",
xlab = "Максимальное количество проданных экземпляров, млн", ylab = "Язык",
main = "Связь максимального количества проданных экземпляров с количеством произведений на языке")
Видно, что с увеличением количества произведений, написанных на каком-либо языке, растёт и максимальные количество проданных экземпляров.
Для того, чтобы визуализировать статистическую характеристику анализируемой совокупности (медиану, нижний и верхний квартили, максимум, минимум, выбросы), используются диаграммы размахов. Например, рассмотрим, как распределено количество проданных экземпляров.
Rython:
plt.figure(figsize=(10, 3))
sns.boxplot(data=df1, x='sales', width=0.4)
plt.xticks(fontsize=9)
plt.xlabel('Количество проданных экземпляров, млн', fontsize=10)
plt.title('Распределение количества проданных экземпляров', fontsize=10, fontweight='bold')
plt.show()
R:
%%R -w 1000 -h 300
boxplot(df1$sales, horizontal = TRUE,
col = "royalblue",
xlab = "Количество проданных экземпляров, млн",
main = "Распределение количества проданных экземпляров")
С помощью диаграмм размахов можно сравнивать распределения разных величин. Для примера рассмотрим, чем отличаются распределения для разных языков.
Python:
tmp2_df = pd.merge(df1, tmp1_df)
tmp2_df = tmp2_df.sort_values(['freq', 'max_sales'], ascending = (True, True))
plt.figure(figsize=(10, 6))
ax = sns.boxplot(data=tmp2_df, x='sales', y='language', dodge=False, hue='freq', palette='deep')
plt.xticks(fontsize=9)
plt.xlabel('Количество проданных экземпляров, млн', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Язык', fontsize=10)
plt.title('Распределение количества проданных экземпляров по языкам', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0,
title='Количество произведений', fontsize=10, title_fontsize=10)
plt.show()
R:
%%R
deep_palette = rep(c("pink", "brown", "purple", "red", "green", "yellow", "royalblue"), times=c(1, 1, 3, 2, 1, 3, 5))
short_deep_palette = c("royalblue", "yellow", "green", "red", "purple", "brown", "pink")
language_levels = c('English', 'Russian', 'French', 'German', 'Japanese', 'Chinese', 'Italian', 'Spanish', 'Hindi', 'Norwegian', 'Swedish', 'Portuguese', 'Dutch', 'Czech', 'Gujarati', 'Yiddish')
%%R
tmp2_df <- data.frame(df1)
tmp2_df$language <- factor(tmp2_df$language, levels = language_levels)
%%R -w 1000 -h 600
# разворачиваем метки осей горизонтально и добавляем дополнительное место справа от графика для легенды и слева для названия оси y
par(las=1, mar=c(5, 7, 4, 16), xpd=TRUE)
boxplot(tmp2_df$sales ~ tmp2_df$language,
horizontal = TRUE, col=deep_palette,
xlab = "Количество проданных экземпляров, млн",
ylab="", # убираем название оси
main = "Распределение количества проданных экземпляров")
# добавляем легенду вне графика
legend("topright", inset=c(-0.25, 0.35), legend=c(1,2,3,4,5,6,132), fill = short_deep_palette, title="Количество произведений")
# устанавливаем назание оси вручную, чтобы оно не накладывалось на метки оси
title(ylab = "Язык", line = 6)
На этих графиках зависимость количества проданных экземпляров от количества книг на языке не настолько явная, как казалось, когда рассматривали только максимумы (см. dotchart).
Для визуализации распределения количественных переменных можно пользоваться одномерными диаграммами рассеивания. Рассмотрим, как распределено количество проданных копий.
Python:
plt.figure(figsize=(10, 3))
sns.stripplot(data=df1, x='sales', orient='h', jitter=True, alpha=0.6, size=7)
plt.xlabel('Количество проданных экземпляров, млн', fontsize=10)
plt.xticks(fontsize=9)
plt.title('Распределение количества проданных экземпляров', fontsize=10, fontweight='bold')
plt.show()
R:
%%R -w 1000 -h 300
stripchart(df1$sales,
xlab = "Количество проданных экземпляров, млн",
vertical = FALSE,
method = "jitter",
pch = 19, col = my_blue, cex = 1.5,
main = "Распределение количества проданных экземпляров")
С помощью диаграмм рассеивания можно сравнивать распределения разных величин. Для примера рассмортим, чем отличаются распределения для разных языков.
Python:
plt.figure(figsize=(10, 6))
ax = sns.stripplot(data=tmp2_df.query("language in ['English', 'Russian', 'French', 'German', 'Japanese', 'Chinese', 'Italian', 'Spanish', 'Hindi', 'Norwegian', 'Swedish', 'Portuguese', 'Dutch', 'Czech', 'Yiddish', 'Gujarati']"),
x='sales', y='language',
orient='h', dodge=False, jitter=True, alpha=0.6, size=7, hue='freq', palette='deep')
plt.xticks(fontsize=9)
plt.xlabel('Количество проданных экземпляров, млн', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Язык', fontsize=10)
plt.title('Распределение количества проданных экземпляров по языкам', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0,
title='Количество произведений', fontsize=10, title_fontsize=10)
plt.show()
R:
%%R -w 1000 -h 600
# разворачиваем метки осей горизонтально и добавляем дополнительное место справа от графика для легенды и слева для названия оси y
par(las=1, mar=c(5, 7, 4, 16), xpd=TRUE)
for(i in 1:length(deep_palette)){
deep_palette[i] <- t_col(deep_palette[i], perc = 60)
}
for(i in 1:length(short_deep_palette)){
short_deep_palette[i] <- t_col(short_deep_palette[i], perc = 60)
}
stripchart(tmp2_df$sales ~ tmp2_df$language,
vertical = FALSE,
method = "jitter",
pch = 19, col = deep_palette, cex = 1.5,
xlab = "Количество проданных экземпляров, млн",
ylab="", # убираем название оси
main = "Распределение количества проданных экземпляров")
# добавляем легенду вне графика
legend("topright", inset=c(-0.25, 0.35), legend=c(1,2,3,4,5,6,132), pch=19, col = short_deep_palette,, title="Количество произведений")
# устанавливаем назание оси вручную, чтобы оно не накладывалось на метки оси
title(ylab = "Язык", line = 6)
На графиках опять наблюдается рост количества проданных экземпляров с ростом количества книг на языке.
Для проверки выбросов с помощью этих критериев формулируются 2 гипотезы (нулевая и альтернативная):
$H_0:$ Подозрительное значение не является выбросом.
$H_1:$ Подозрительное значение - выброс.
Если полученый $pvalue$ меньше заданного заранее уровня значимости $\alpha$ (возьмём $\alpha = 0.05$), то $H_0$ должна быть отвергнута в пользу $H_1$, иначе - нет оснований отвергнуть нулевую гипотезу.
Также и критерий Граббса, и Q-тест Диксона основаны на предположении о нормальности распределения. Поэтому применим их к данным, близким к нормальным (см. график ядерных оценок плотности) - распределению годов первой публикации.
Python:
В реализации теста Граббса на Python $pvalue$ "спрятан", поэтому выводим только результат.
min_out_index = grubbs.min_test_indices(df1.first_published, alpha=0.05)
max_out_index = grubbs.max_test_indices(df1.first_published, alpha=0.05)
df1.iloc[min_out_index + max_out_index]
| title | author | gender | birth_year | age | language | first_published | sales | genres | |
|---|---|---|---|---|---|---|---|---|---|
| 31 | The Divine Comedy (La Divina Commedia) | Dante Alighieri | M | 1265 | 39 | Italian | 1304 | 11.5 | [Classics, Poetry, Fiction, Literature, Philos... |
R:
%%R
res <- grubbs.test(df1[]$first_published, type = 10)
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет")
Grubbs test for one outlier data: df1[]$first_published G = 10.27459, U = 0.38981, p-value < 2.2e-16 alternative hypothesis: lowest value 1304 is an outlier Есть ли основания отбросить гипотезу? Да
%%R
res <- grubbs.test(df1[]$first_published, type = 10, opposite = TRUE)
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет")
Grubbs test for one outlier data: df1[]$first_published G = 0.86350, U = 0.99569, p-value = 1 alternative hypothesis: highest value 2018 is an outlier Есть ли основания отбросить гипотезу? Нет
%%R
df1[df1$first_published == 1304, ]
title author gender birth_year age
32 The Divine Comedy (La Divina Commedia) Dante Alighieri M 1265 39
language first_published sales
32 Italian 1304 11.5
genres
32 Classics, Poetry, Fiction, Literature, Philosophy, Religion, Fantasy
Тест Граббса нашёл один выброс - "Божественную комедию" Данте Алигьери (1304).
Тест Диксона может обрабатывать выборки размером не более 30. Поэтому выделим диапазон из не более чем 30 значений, в который потенциально может попасть выброс. Также если в диапазон не попала часть одинаковых значений, удалим все такие значения.
Python:
min_df = df1.sort_values(['first_published']).reset_index()[0:31]
min_df = min_df[min_df.first_published != min_df.first_published[30]]
print('Количество элементов для проверки минимума:', min_df.shape[0])
Количество элементов для проверки минимума: 30
max_df = df1.sort_values(['first_published'], ascending=False).reset_index()[0:31]
max_df = max_df[max_df.first_published != max_df.first_published[30]]
print('Количество элементов для проверки максимума:', max_df.shape[0])
Количество элементов для проверки максимума: 29
R:
%%R
min_df <- df1[order(df1$first_published, decreasing=FALSE), ]$first_published[1:31]
min_df <- min_df[min_df != min_df[31]]
cat("Количество элементов для проверки минимума:", length(min_df))
Количество элементов для проверки минимума: 30
%%R
max_df <- df1[order(df1$first_published, decreasing=TRUE), ]$first_published[1:31]
max_df <- max_df[max_df != max_df[31]]
cat("Количество элементов для проверки максимума:", length(max_df))
Количество элементов для проверки максимума: 29
Для минимальных и максимальных выбросов получились выборки размером 30 и 29 соответственно.
Python:
На Python нет готовой реализации теста Диксона, поэтому сделаем её сами. Для выборок размера 30 и 29 при уровне значимости 0.95 критические значения равны 0.290 и 0.301 соответственно (это табличные данные, полученные моделированием; численно получить их можно только для выборки размера $n=3$).
out = []
maX = min_df.first_published[29]
miN = min_df.first_published[0]
second_min = min(min_df[min_df.first_published != miN].first_published)
min_q = (second_min - miN) / (maX - miN)
if (min_q >= 0.290):
out.append(miN)
maX = max_df.first_published[0]
second_max = max(max_df[max_df.first_published != maX].first_published)
miN = max_df.first_published[28]
max_q = (maX - second_max) / (maX - miN)
if (max_q >= 0.310):
out.append(maX)
df1[df1.first_published.isin(out)]
| title | author | gender | birth_year | age | language | first_published | sales | genres | |
|---|---|---|---|---|---|---|---|---|---|
| 31 | The Divine Comedy (La Divina Commedia) | Dante Alighieri | M | 1265 | 39 | Italian | 1304 | 11.5 | [Classics, Poetry, Fiction, Literature, Philos... |
R:
%%R
res <- dixon.test(min_df)
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет")
Dixon test for outliers data: min_df Q = 0.77057, p-value < 2.2e-16 alternative hypothesis: lowest value 1304 is an outlier Есть ли основания отбросить гипотезу? Да
%%R
res <- dixon.test(max_df)
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет")
Dixon test for outliers data: max_df Q = 0.13333, p-value = 0.6394 alternative hypothesis: highest value 2018 is an outlier Есть ли основания отбросить гипотезу? Нет
Тест Диксона дал тот же результат, что и тест Граббса. Визуализируем его.
Python:
plt.figure(figsize=(10, 3))
ax = sns.stripplot(data=df1[df1.first_published.isin(out)], x='first_published', orient='h', jitter=True,
alpha=0.6, size=7, color='red', label='Выброс')
sns.stripplot(data=df1[~df1.first_published.isin(out)], x='first_published', orient='h', jitter=True,
alpha=0.6, size=7, color='green', label='Остальная выборка')
plt.xticks(fontsize=9)
plt.xlabel('Год первой публикации', fontsize=10)
plt.title('Распределение первых публикаций', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0,
title='')
plt.show()
R:
%%R -w 1000 -h 300
my_green <- t_col("forestgreen", perc = 60)
stripchart(df1[df1$first_published != 1304, ]$first_published,
method = "jitter",
vertical = FALSE,
pch = 19, col = my_green, cex = 1.5,
xlim = c(1304, 2018),
xlab = "Распределение первых публикаций",
main = "Распределение количества проданных экземпляров")
stripchart(df1[df1$first_published == 1304, ]$first_published,
method = "jitter",
pch = 19, col = my_red, cex = 1.5, add = TRUE)
Рассмотрим распределение возрастов авторов.
Python:
plt.figure(figsize=(10, 3))
sns.stripplot(data=df1, x='age', orient='h', jitter=True,
alpha=0.6, size=7)
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.title('Распределение возрастов авторов', fontsize=10, fontweight='bold')
plt.show()
R:
%%R -w 1000 -h 300
stripchart(df1$age,
xlab = "Возраст",
vertical = FALSE,
method = "jitter",
pch = 19, col = my_blue, cex = 1.5,
main = "Распределение возрастов авторов")
Удалим возраст в некоторых строках.
Python:
df_copy = df1.sort_values(['age', 'author']).reset_index()
df_copy.iloc[[36, 89, 141, 147]]
| index | title | author | gender | birth_year | age | language | first_published | sales | genres | |
|---|---|---|---|---|---|---|---|---|---|---|
| 36 | 104 | Anne of Green Gables | Lucy Maud Montgomery | F | 1874 | 34 | English | 1908 | 50.0 | [Classics, Fiction, Young Adult, Historical Fi... |
| 89 | 165 | Man's Search for Meaning (Ein Psychologe erleb... | Viktor Frankl | M | 1905 | 41 | German | 1946 | 12.0 | [Nonfiction, Psychology, Philosophy, History, ... |
| 141 | 18 | The Lion, the Witch and the Wardrobe | C. S. Lewis | M | 1898 | 52 | English | 1950 | 85.0 | [Fantasy, Classics, Fiction, Young Adult, Chil... |
| 147 | 85 | Heidi | Johanna Spyri | F | 1827 | 53 | German | 1880 | 50.0 | [Classics, Fiction, Childrens, Young Adult, Hi... |
df_copy['age'][36] = np.nan
df_copy['age'][89] = np.nan
df_copy['age'][141] = np.nan
df_copy['age'][147] = np.nan
df_copy['age'].isna().sum()
4
R:
%%R
df_copy <- df1[order(c(df1$age)), ]
df_copy[c(37, 90, 142, 148), ]
title
105 Anne of Green Gables
166 Man's Search for Meaning (Ein Psychologe erlebt das Konzentrationslager)
19 The Lion, the Witch and the Wardrobe
86 Heidi
author gender birth_year age language first_published sales
105 Lucy Maud Montgomery F 1874 34 English 1908 50
166 Viktor Frankl M 1905 41 German 1946 12
19 C. S. Lewis M 1898 52 English 1950 85
86 Johanna Spyri F 1827 53 German 1880 50
genres
105 Classics, Fiction, YoungAdult, HistoricalFiction, Childrens, MiddleGrade, Audiobook
166 Nonfiction, Psychology, Philosophy, History, SelfHelp, Memoir, Biography
19 Fantasy, Classics, Fiction, YoungAdult, Childrens, MiddleGrade, Adventure
86 Classics, Fiction, Childrens, YoungAdult, HistoricalFiction, MiddleGrade, Historical
%%R
df_copy$age[37] <- NA
df_copy$age[90] <- NA
df_copy$age[142] <- NA
df_copy$age[148] <- NA
sum(is.na(df_copy$age))
[1] 4
Заполним пропуски средним по всем значениям и cравним получившееся распределение с истинными.
Python:
# одинаковые границы отрезков разбиения для гистограмм
edges = np.histogram_bin_edges(df1['age'], bins=14, range=(15,85))
edges
array([15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75.,
80., 85.])
all_mean_df = df_copy.copy()
all_mean = all_mean_df.mean()['age']
all_mean_df['age'][36] = all_mean
all_mean_df['age'][89] = all_mean
all_mean_df['age'][141] = all_mean
all_mean_df['age'][147] = all_mean
print(all_mean)
42.7953216374269
fig = plt.figure(figsize=(10,3))
ax1 = plt.subplot(121)
sns.histplot(data=df1.age, kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax1.set_title('Истинные значения', fontsize=10, fontweight='bold')
ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
sns.histplot(data=all_mean_df.age, kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax2.set_title('Пропуски заполнены средним по всем значениям', fontsize=10, fontweight='bold')
plt.show()
R:
%%R
# одинаковые границы отрезков разбиения для гистограмм
edges = hist(df1$age, freq = FALSE, breaks = "FD", plot = FALSE)
edges[1]
$breaks [1] 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85
%%R
all_mean_df <- data.frame(df_copy)
all_mean <- mean(all_mean_df$age, na.rm = TRUE)
all_mean_df$age[37] <- all_mean
all_mean_df$age[90] <- all_mean
all_mean_df$age[142] <- all_mean
all_mean_df$age[148] <- all_mean
cat(all_mean)
42.79532
%%R -w 1000 -h 300
par(mfrow = c(1,2))
hist(df1$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "Истинные значения")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 3)
hist(all_mean_df$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "Пропуски заполнены средним по всем значениям")
lines(density(all_mean_df$age, na.rm = TRUE, adjust=), col = "royalblue", lwd = 3)
Заполним пропуски медианой по всем значениям и вычислим отклонение заполнений от истинных значений.
Python:
all_median_df = df_copy.copy()
all_median = all_median_df.median()['age']
all_median_df['age'][36] = all_median
all_median_df['age'][89] = all_median
all_median_df['age'][141] = all_median
all_median_df['age'][147] = all_median
print(all_median)
41.0
fig = plt.figure(figsize=(10,3))
ax1 = plt.subplot(121)
sns.histplot(data=df1.age, kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax1.set_title('Истинные значения', fontsize=10, fontweight='bold')
ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
sns.histplot(data=all_median_df.age, kde=True, bins=edges, alpha=0.3, stat='density')
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax2.set_title('Пропуски заполнены медианой по всем значениям', fontsize=10, fontweight='bold')
plt.show()
R:
%%R
all_median_df <- data.frame(df_copy)
all_median <- median(all_median_df$age, na.rm = TRUE)
all_median_df$age[37] <- all_median
all_median_df$age[90] <- all_median
all_median_df$age[142] <- all_median
all_median_df$age[148] <- all_median
cat(all_median)
41
%%R -w 1000 -h 300
par(mfrow = c(1,2))
hist(df1$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "Истинные значения")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 3)
hist(all_median_df$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "Пропуски заполнены медианой по всем значениям")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "royalblue", lwd = 3)
Пропуски можно заполнять средними и медианами не по всем значениям, а по категории. Для примера рассмотрим, чем отличаются распределения для женщин и мужчин.
Python:
plt.figure(figsize=(10, 3))
ax = sns.stripplot(data=df1, x='age', orient='h', jitter=True, dodge=True, hue='gender',
alpha=0.4, size=7, palette=['red', 'blue'])
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.title('Распределение возрастов авторов', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0,
title='Пол автора', fontsize=10, title_fontsize=10)
plt.show()
R:
%%R -w 1000 -h 300
# добавляем дополнительное место справа от графика для легенды
par(mar=c(5, 4, 4, 16), xpd=TRUE)
stripchart(df1$age ~ df1$gender,
vertical = FALSE,
method = "jitter",
pch = 19, col = c(my_red, my_blue), cex = 1.5,
xlab = "Возраст",
ylab="", yaxt="n",
main = "Распределение возрастов авторов")
# добавляем легенду вне графика
legend("topright", inset=c(-0.15, 0.35), legend=c("M", "F"), pch=19,
col = c(my_blue, my_red), title = "Пол автора")
Заполним пропуски средними по категориям и вычислим отклонение заполнений от истинных значений.
Python:
gender_mean_df = df_copy.copy()
m_mean = gender_mean_df[gender_mean_df.gender == 'M'].mean()['age']
f_mean = gender_mean_df[gender_mean_df.gender == 'F'].mean()['age']
gender_mean_df['age'][36] = m_mean if (gender_mean_df['gender'][36] == 'M') else f_mean
gender_mean_df['age'][89] = m_mean if (gender_mean_df['gender'][89] == 'M') else f_mean
gender_mean_df['age'][141] = m_mean if (gender_mean_df['gender'][141] == 'M') else f_mean
gender_mean_df['age'][147] = m_mean if (gender_mean_df['gender'][147] == 'M') else f_mean
print(m_mean, f_mean)
43.4 41.55357142857143
fig = plt.figure(figsize=(10,9))
ax1 = plt.subplot(321)
sns.histplot(data=df1.age,
kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax1.set_title('Истинные значения', fontsize=10, fontweight='bold')
ax2 = plt.subplot(322, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_mean_df.age,
kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax2.set_title('Пропуски заполнены средними по категориям', fontsize=10, fontweight='bold')
ax3 = plt.subplot(323, sharex=ax1, sharey=ax1)
sns.histplot(data=df1[df1.gender == 'M'].age,
kde=True, bins=edges, alpha=0.3, color='blue', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax3.set_title('Истинные значения (муж)', fontsize=10, fontweight='bold')
ax4 = plt.subplot(324, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_mean_df[gender_mean_df.gender == 'M'].age,
kde=True, bins=edges, alpha=0.3, color='blue', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax4.set_title('Пропуски заполнены средними по категориям (муж)', fontsize=10, fontweight='bold')
ax5 = plt.subplot(325, sharex=ax1, sharey=ax1)
sns.histplot(data=df1.age[df1.gender == 'F'],
kde=True, bins=edges, alpha=0.3, color='red', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax5.set_title('Истинные значения (жен)', fontsize=10, fontweight='bold')
ax6 = plt.subplot(326, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_mean_df[gender_mean_df.gender == 'F'].age,
kde=True, bins=edges, alpha=0.3, color='red', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax6.set_title('Пропуски заполнены средними по категориям (жен)', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=2)
plt.show()
R:
%%R
gender_mean_df <- data.frame(df_copy)
m_mean <- mean(gender_mean_df[gender_mean_df$gender %like% "M", "age"], na.rm = TRUE)
f_mean <- mean(gender_mean_df[gender_mean_df$gender %like% "F", "age"], na.rm = TRUE)
gender_mean_df$age[37] <- if (gender_mean_df$gender[37] == "M") m_mean else f_mean
gender_mean_df$age[90] <- if (gender_mean_df$gender[90] == "M") m_mean else f_mean
gender_mean_df$age[142] <- if (gender_mean_df$gender[142] == "M") m_mean else f_mean
gender_mean_df$age[148] <- if (gender_mean_df$gender[148] == "M") m_mean else f_mean
cat(m_mean, f_mean)
43.4 41.55357
%%R -w 1000 -h 900
par(mfrow = c(3,2))
hist(df1$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "Истинные значения")
lines(density(df1$age, adjust=), col = "blue", lwd = 3)
hist(gender_mean_df$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "Пропуски заполнены средними по категориям")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "blue", lwd = 3)
hist(df1[df1$gender %like% "M", ]$age, freq = FALSE, breaks = edges$breaks,
col = my_blue,
xlab = "Возраст", ylab = "",
main = "Истинные значения (муж)")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 3)
hist(gender_mean_df[gender_mean_df$gender %like% "M", ]$age, freq = FALSE, breaks = edges$breaks,
col = my_blue,
xlab = "Возраст", ylab = "",
main = "Пропуски заполнены средними по категориям (муж)")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "royalblue", lwd = 3)
hist(df1[df1$gender %like% "F", ]$age, freq = FALSE, breaks = edges$breaks,
col = my_red,
xlab = "Возраст", ylab = "",
main = "Истинные значения (жен)")
lines(density(df1$age, adjust=), col = "red", lwd = 3)
hist(gender_mean_df[gender_mean_df$gender %like% "F", ]$age, freq = FALSE, breaks = edges$breaks,
col = my_red,
xlab = "Возраст", ylab = "",
main = "Пропуски заполнены средними по категориям (жен)")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "red", lwd = 3)
Заполним пропуски медианами по категориям и вычислим отклонение заполнений от истинных значений.
Python:
gender_median_df = df_copy.copy()
m_median = gender_median_df[gender_median_df.gender == 'M'].median()['age']
f_median = gender_median_df[gender_median_df.gender == 'F'].median()['age']
gender_median_df['age'][36] = m_median if (gender_median_df['gender'][36] == 'M') else f_median
gender_median_df['age'][89] = m_median if (gender_median_df['gender'][89] == 'M') else f_median
gender_median_df['age'][141] = m_median if (gender_median_df['gender'][141] == 'M') else f_median
gender_median_df['age'][147] = m_median if (gender_median_df['gender'][147] == 'M') else f_median
print(m_median, f_median)
41.0 40.0
fig = plt.figure(figsize=(10,9))
ax1 = plt.subplot(321)
sns.histplot(data=df1.age,
kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax1.set_title('Истинные значения', fontsize=10, fontweight='bold')
ax2 = plt.subplot(322, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_median_df.age,
kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax2.set_title('Пропуски заполнены медианами по категориям', fontsize=10, fontweight='bold')
ax3 = plt.subplot(323, sharex=ax1, sharey=ax1)
sns.histplot(data=df1[df1.gender == 'M'].age,
kde=True, bins=edges, alpha=0.3, color='blue', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax3.set_title('Истинные значения (муж)', fontsize=10, fontweight='bold')
ax4 = plt.subplot(324, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_median_df[gender_median_df.gender == 'M'].age,
kde=True, bins=edges, alpha=0.3, color='blue', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax4.set_title('Пропуски заполнены медианами по категориям (муж)', fontsize=10, fontweight='bold')
ax5 = plt.subplot(325, sharex=ax1, sharey=ax1)
sns.histplot(data=df1.age[df1.gender == 'F'],
kde=True, bins=edges, alpha=0.3, color='red', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax5.set_title('Истинные значения (жен)', fontsize=10, fontweight='bold')
ax6 = plt.subplot(326, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_median_df[gender_median_df.gender == 'F'].age,
kde=True, bins=edges, alpha=0.3, color='red', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax6.set_title('Пропуски заполнены медианами по категориям (жен)', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=2)
plt.show()
R:
%%R
gender_median_df <- data.frame(df_copy)
m_median <- median(gender_median_df[gender_median_df$gender %like% "M", "age"], na.rm = TRUE)
f_median <- median(gender_median_df[gender_median_df$gender %like% "F", "age"], na.rm = TRUE)
gender_median_df$age[37] <- if (gender_median_df$gender[37] == "M") m_median else f_median
gender_median_df$age[90] <- if (gender_median_df$gender[90] == "M") m_median else f_median
gender_median_df$age[142] <- if (gender_median_df$gender[142] == "M") m_median else f_median
gender_median_df$age[148] <- if (gender_median_df$gender[148] == "M") m_median else f_median
cat(m_median, f_median)
41 40
%%R -w 1000 -h 900
par(mfrow = c(3,2))
hist(df1$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "Истинные значения")
lines(density(df1$age, adjust=), col = "blue", lwd = 3)
hist(gender_median_df$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "Пропуски заполнены медианами по категориям")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "blue", lwd = 3)
hist(df1[df1$gender %like% "M", ]$age, freq = FALSE, breaks = edges$breaks,
col = my_blue,
xlab = "Возраст", ylab = "",
main = "Истинные значения (муж)")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 3)
hist(gender_median_df[gender_median_df$gender %like% "M", ]$age, freq = FALSE, breaks = edges$breaks,
col = my_blue,
xlab = "Возраст", ylab = "",
main = "Пропуски заполнены медианами по категориям (муж)")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "royalblue", lwd = 3)
hist(df1[df1$gender %like% "F", ]$age, freq = FALSE, breaks = edges$breaks,
col = my_red,
xlab = "Возраст", ylab = "",
main = "Истинные значения (жен)")
lines(density(df1$age, adjust=), col = "red", lwd = 3)
hist(gender_median_df[gender_median_df$gender %like% "F", ]$age, freq = FALSE, breaks = edges$breaks,
col = my_red,
xlab = "Возраст", ylab = "",
main = "Пропуски заполнены медианами по категориям (жен)")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "red", lwd = 3)
Сравним разные способы заполнения.
Python:
real_age = df1.sort_values(['age']).reset_index()[['age']].iloc[[36, 89, 141, 147]]
all_mean_age = all_mean_df.reset_index()[['age']].iloc[[36, 89, 141, 147]]
all_median_age = all_median_df.reset_index()[['age']].iloc[[36, 89, 141, 147]]
gender_mean_age = gender_mean_df.reset_index()[['age']].iloc[[36, 89, 141, 147]]
gender_median_age = gender_median_df.reset_index()[['age']].iloc[[36, 89, 141, 147]]
all_mean_age['age'] = real_age['age'] - all_mean_age['age']
all_median_age['age'] = real_age['age'] - all_median_age['age']
gender_mean_age['age'] = real_age['age'] - gender_mean_age['age']
gender_median_age['age'] = real_age['age'] - gender_median_age['age']
fig = plt.figure(figsize=(10,6))
ax1 = plt.subplot(221)
plt.bar(data=all_mean_age.query("index in ['36', '89', '141', '147']"), x=['36', '89', '141', '147'],
height=all_mean_age.age)
plt.xticks(fontsize=9)
plt.xlabel('Индекс пропуска', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Разница с истинным значением', fontsize=10)
plt.title('Пропуски заполнены средним по всем', fontsize=10, fontweight='bold')
ax2 = plt.subplot(222, sharex=ax1, sharey=ax1)
plt.bar(data=all_median_age.query("index in ['36', '89', '141', '147']"), x=['36', '89', '141', '147'],
height=all_median_age.age)
plt.xticks(fontsize=9)
plt.xlabel('Индекс пропуска', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Разница с истинным значением', fontsize=10)
plt.title('Пропуски заполнены медианой по всем', fontsize=10, fontweight='bold')
ax3 = plt.subplot(223, sharex=ax1, sharey=ax1)
plt.bar(data=gender_mean_age.query("index in ['36', '89', '141', '147']"), x=['36', '89', '141', '147'],
height=gender_mean_age.age)
plt.xticks(fontsize=9)
plt.xlabel('Индекс пропуска', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Разница с истинным значением', fontsize=10)
plt.title('Пропуски заполнены средними по категориям', fontsize=10, fontweight='bold')
ax4 = plt.subplot(224, sharex=ax1, sharey=ax1)
plt.bar(data=gender_median_age.query("index in ['36', '89', '141', '147']"), x=['36', '89', '141', '147'],
height=gender_median_age.age)
plt.xticks(fontsize=9)
plt.xlabel('Индекс пропуска', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Разница с истинным значением', fontsize=10)
plt.title('Пропуски заполнены медианами по категориям', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=3)
plt.show()
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
R:
%%R
real_age <- df1[order(c(df1$age)), ]$age[c(37, 90, 142, 148)]
all_mean_age <- all_mean_df$age[c(37, 90, 142, 148)]
all_median_age <- all_median_df$age[c(37, 90, 142, 148)]
gender_mean_age <- gender_mean_df$age[c(37, 90, 142, 148)]
gender_median_age <- gender_median_df$age[c(37, 90, 142, 148)]
all_mean_age <- real_age - all_mean_age
all_median_age<- real_age - all_median_age
gender_mean_age <- real_age - gender_mean_age
gender_median_age <- real_age - gender_median_age
%%R -w 1000 -h 600
par(mfrow = c(2,2))
barplot(all_mean_age, names.arg = c("37", "90", "142", "148"),
col = "lightblue",
xlab = "Индекс пропуска", ylab = "Разница с истинным значением",
main = "Пропуски заполнены средним по всем")
barplot(all_median_age, names.arg = c("37", "90", "142", "148"),
col = "lightblue",
xlab = "Индекс пропуска", ylab = "Разница с истинным значением",
main = "Пропуски заполнены медианой по всем")
barplot(gender_mean_age, names.arg = c("37", "90", "142", "148"),
col = "lightblue",
xlab = "Индекс пропуска", ylab = "Разница с истинным значением",
main = "Пропуски заполнены средними по категориям")
barplot(gender_median_age, names.arg = c("37", "90", "142", "148"),
col = "lightblue",
xlab = "Индекс пропуска", ylab = "Разница с истинным значением",
main = "Пропуски заполнены медианами по категориям")
В данном случае разные способы заполнения пропусков дают примерно одинаковый результат, но при заполнениях медианами значения получаются меньше из-за ассиметрии данных.
Сгенерируем данные из нормальных распределений с различными параметрами:
Python:
data1 = np.random.normal(0, 1, size=100)
data2 = np.random.normal(0, 2, size=100)
data3 = np.random.normal(0, 1, size=1000)
data4 = np.random.normal(2, 1, size=1000)
R:
%%R
data1 <- rnorm(100, mean = 0, sd = 1)
data2 <- rnorm(100, mean = 0, sd = 2)
data3 <- rnorm(1000, mean = 0, sd = 1)
data4 <- rnorm(1000, mean = 2, sd = 1)
Для проверки нормальности можно использовать эмпирические функции распределения, сравнивая их графики с графиками теоретических функций распределения. В случае нормального распределения они должны быть похожи.
Python:
fig = plt.figure(figsize=(10,6))
x = np.linspace(-10, 10, 1000)
ax1 = plt.subplot(221)
plt.plot(x, stats.norm.cdf(x), color='red')
sns.ecdfplot(data=data1)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax1.set_title('N(0,1), 100 наблюдений', fontsize=10, fontweight='bold')
ax2 = plt.subplot(222, sharex=ax1, sharey=ax1)
plt.plot(x, stats.norm.cdf(x, 0, 2), color='red')
sns.ecdfplot(data=data2)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax2.set_title('N(0,2), 100 наблюдений', fontsize=10, fontweight='bold')
ax3 = plt.subplot(223, sharex=ax1, sharey=ax1)
plt.plot(x, stats.norm.cdf(x), color='red')
sns.ecdfplot(data=data3)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax3.set_title('N(0,1), 1000 наблюдений', fontsize=10, fontweight='bold')
ax4 = plt.subplot(224, sharex=ax1, sharey=ax1)
plt.plot(x, stats.norm.cdf(x, 2, 1), color='red')
sns.ecdfplot(data=data4)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax4.set_title('N(2,1), 1000 наблюдений', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=2)
plt.show()
R:
%%R -w 1000 -h 600
par(mfrow = c(2,2))
plot(sort(data1), pnorm(sort(data1), mean = 0, sd = 1), type = "l",
xlim = c(-10, 10), ylim = c(-0.01, 1.01),
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="N(0,1), 100 наблюдений")
plot(ecdf(data1), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
plot(sort(data2), pnorm(sort(data2), mean = 0, sd = 2), type = "l",
xlim = c(-10, 10), ylim = c(-0.01, 1.01),
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="N(0,2), 100 наблюдений")
plot(ecdf(data2), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
plot(sort(data3), pnorm(sort(data3), mean = 0, sd = 1), type = "l",
xlim = c(-10, 10), ylim = c(-0.01, 1.01),
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="N(0,1), 1000 наблюдений")
plot(ecdf(data3), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
plot(sort(data4), pnorm(sort(data4), mean = 2, sd = 1), type = "l",
xlim = c(-10, 10), ylim = c(-0.01, 1.01),
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="N(2,1), 1000 наблюдений")
plot(ecdf(data4), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
Проверять, получены ли данные из некоторого теоретического распределения, можно с помощью графиков QQ ("квантиль-квантиль"). В случае нормальных распределений с различными параметрами точки должны лежать вдоль соответствующих прямых.
Python:
fig = plt.figure(figsize=(10,6))
ax1 = plt.subplot(221)
sm.qqplot(data1, ax=ax1, line='45')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
ax1.set_xlim(-6, 6)
plt.yticks(fontsize=9)
ax1.set_ylim(-6, 6)
plt.ylabel('Выборочные квантили', fontsize=10)
ax1.set_title('N(0,1), 100 наблюдений', fontsize=10, fontweight='bold')
ax2 = plt.subplot(222, sharex=ax1, sharey=ax1)
sm.qqplot(data2, ax=ax2, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax2.set_title('N(0,2), 100 наблюдений', fontsize=10, fontweight='bold')
ax3 = plt.subplot(223, sharex=ax1, sharey=ax1)
sm.qqplot(data3, ax=ax3, line='45')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax3.set_title('N(0,1), 1000 наблюдений', fontsize=10, fontweight='bold')
ax4 = plt.subplot(224, sharex=ax1, sharey=ax1)
sm.qqplot(data4, ax=ax4, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax4.set_title('N(2,1), 1000 наблюдений', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=2)
plt.show()
R:
%%R -w 1000 -h 600
par(mfrow = c(2,2))
qqnorm(data1,
xlim = c(-6, 6), ylim = c(-6, 6),
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "N(0,1), 100 наблюдений")
qqline(data1,
col = "red", lwd = 2)
qqnorm(data2,
xlim = c(-6, 6), ylim = c(-6, 6),
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "N(0,2), 100 наблюдений")
qqline(data2,
col = "red", lwd = 2)
qqnorm(data3,
xlim = c(-6, 6), ylim = c(-6, 6),
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "N(0,1), 1000 наблюдений")
qqline(data3,
col = "red", lwd = 2)
qqnorm(data4,
xlim = c(-6, 6), ylim = c(-6, 6),
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "N(2,1), 1000 наблюдений")
qqline(data4,
col = "red", lwd = 2)
Ещё один способ визуально проверить нормальность распределения - метод огибающих.
Python:
def envelmet(x, title, ax):
# точки
z = (x - x.mean()) / x.std()
x_qq, _ = stats.probplot(z)
ax.plot(x_qq[0], x_qq[1], '.k')
# линия
ax.plot((-3, 3), (-3, 3), c='red')
ax.tick_params(axis='both', labelsize=9)
ax.set_xlabel('Квантили', fontsize=10)
ax.set_ylabel('Статистики', fontsize=10)
ax.set_title(title, fontsize=10, fontweight='bold')
fig, ax = plt.subplots(nrows=2, ncols=2, figsize = (10, 6))
envelmet(data1, 'N(0,1), 100 наблюдений', ax[0, 0])
envelmet(data2, 'N(0,2), 100 наблюдений', ax[0, 1])
envelmet(data3, 'N(0,1), 1000 наблюдений', ax[1, 0])
envelmet(data4, 'N(2,1), 1000 наблюдений', ax[1, 1])
fig.tight_layout()
plt.show()
R:
%%R -w 1000 -h 600
par(mfrow = c(2,2))
envelmet <- function(x, title){
# точки
z <- (x - mean(x)) / sqrt(var(x))
x_qq <- qqnorm(z, plot.it = FALSE)
x_qq <- lapply(x_qq, sort)
plot(x_qq,
xlim = c(-3, 3), ylim = c(-3, 3),
xlab = "Квантили", ylab = "Статистики",
main = title)
# огибающие
x.gen <- function(dat, mle) rnorm(length(dat))
x.qqboot <- boot(z, sort, R = 999,
sim = "parametric", ran.gen = x.gen)
sapply(1:999, function(i) lines(x_qq$x, x.qqboot$t[i,], type = "l", col = "grey"))
points (x_qq, pch = 20)
#прямая
lines(c(-3, 3), c(-3, 3), col = "red", lwd = 2)
}
envelmet(data1, "N(0,1), 100 наблюдений")
envelmet(data2, "N(0,2), 100 наблюдений")
envelmet(data3, "N(0,1), 1000 наблюдений")
envelmet(data4, "N(2,1), 1000 наблюдений")
Для каждого из критериев формулируется нулевая гипотеза $H_0$ о характере распределения выборки $X_1, \dots, X_n$. Если полученный при использовании критерия $pvalue$ меньше заданного заранее уровня значимости $\alpha$ (положим $\alpha = 0.5$), то $H_0$ должна быть отвергнута, иначе - нет оснований отбросить нулевую гипотезу.
т.е.данные получены из распределения с заданными параметрами.
Положим $F(x, \theta_0) \equiv \Phi(x)$, чтобы проверить, является ли распределение стандартным нормальным.
Python:
_, pvalue = kstest(data1, 'norm')
print('N(0,1), 100 наблюдений:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < 0.05) else 'Нет')
_, pvalue = kstest(data2, 'norm')
print('N(0,2), 100 наблюдений:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < 0.05) else 'Нет')
_, pvalue = kstest(data3, 'norm')
print('N(0,1), 1000 наблюдений:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < 0.05) else 'Нет')
_, pvalue = kstest(data4, 'norm')
print('N(2,1), 1000 наблюдений:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
pvalue = 0.48487951108184757
Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
pvalue = 0.0002731978404077222
Есть ли основания отбросить гипотезу? Да
N(0,1), 1000 наблюдений:
pvalue = 0.6034541166986684
Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
pvalue = 0.0
Есть ли основания отбросить гипотезу? Да
R:
%%R
res <- ks.test(data1, "pnorm")
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- ks.test(data2, "pnorm")
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- ks.test(data3, "pnorm")
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- ks.test(data4, "pnorm")
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------ Asymptotic one-sample Kolmogorov-Smirnov test data: data1 D = 0.071778, p-value = 0.6815 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Нет ------------------ N(0,2), 100 наблюдений ------------------ Asymptotic one-sample Kolmogorov-Smirnov test data: data2 D = 0.22247, p-value = 0.0001005 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Да ------------------ N(0,1), 1000 наблюдений ------------------ Asymptotic one-sample Kolmogorov-Smirnov test data: data3 D = 0.014548, p-value = 0.984 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Нет ------------------ N(2,1), 1000 наблюдений ------------------ Asymptotic one-sample Kolmogorov-Smirnov test data: data4 D = 0.69077, p-value < 2.2e-16 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Да
т.е. выборка имеет нормальное распределение, параметры которого заранее не известны.
Также для критерия Шапиро-Уилка необходимо, чтобы размер выборки был не более 5000.
Python:
_, pvalue = shapiro(data1)
print('N(0,1), 100 наблюдений:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
_, pvalue = shapiro(data2)
print('N(0,2), 100 наблюдений:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
_, pvalue = shapiro(data3)
print('N(0,1), 1000 наблюдений:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
_, pvalue = shapiro(data4)
print('N(2,1), 1000 наблюдений:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
pvalue = 0.6277329921722412
Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
pvalue = 0.04325342923402786
Есть ли основания отбросить гипотезу? Да
N(0,1), 1000 наблюдений:
pvalue = 0.9437151551246643
Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
pvalue = 0.2394869327545166
Есть ли основания отбросить гипотезу? Нет
R:
%%R
res <- shapiro.test(data1)
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- shapiro.test(data2)
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- shapiro.test(data3)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- shapiro.test(data4)
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------ Shapiro-Wilk normality test data: data1 W = 0.99228, p-value = 0.841 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,2), 100 наблюдений ------------------ Shapiro-Wilk normality test data: data2 W = 0.99281, p-value = 0.8764 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,1), 1000 наблюдений ------------------ Shapiro-Wilk normality test data: data3 W = 0.99847, p-value = 0.5338 Есть ли основания отбросить гипотезу? Нет ------------------ N(2,1), 1000 наблюдений ------------------ Shapiro-Wilk normality test data: data4 W = 0.99719, p-value = 0.0786 Есть ли основания отбросить гипотезу? Нет
т.е. выборка имеет нормальное распределение, параметры которого заранее не известны.
Python:
Так как в Python реализация критерия Андерсона-Дарлинга (scipy.stats.anderson()) не возвращает pvalue, нужно сравнить полученную тестовую статистику с табличным критическим значением (для уровня значимости $\alpha = 0.05$ и выборок размера 100 и 1000 критическое значение равно 0.759 и 0.784 соответственно).
res = anderson(data1, 'norm')
print('N(0,1), 100 наблюдений')
print(' statistic =', res.statistic)
print(' Есть ли основания отбросить гипотезу?',
'Да' if (res.critical_values[np.where(res.significance_level == 5.)] < res.statistic) else 'Нет')
res = anderson(data2, 'norm')
print('N(0,2), 100 наблюдений')
print(' statistic =', res.statistic)
print(' Есть ли основания отбросить гипотезу?',
'Да' if (res.critical_values[np.where(res.significance_level == 5.)] < res.statistic) else 'Нет')
res = anderson(data3, 'norm')
print('N(0,1), 1000 наблюдений')
print(' statistic =', res.statistic)
print(' Есть ли основания отбросить гипотезу?',
'Да' if (res.critical_values[np.where(res.significance_level == 5.)] < res.statistic) else 'Нет')
res = anderson(data4, 'norm')
print('N(2,1), 1000 наблюдений')
print(' statistic =', res.statistic)
print(' Есть ли основания отбросить гипотезу?',
'Да' if (res.critical_values[np.where(res.significance_level == 5.)] < res.statistic) else 'Нет')
N(0,1), 100 наблюдений
statistic = 0.3803493043023565
Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений
statistic = 0.5864913565249026
Есть ли основания отбросить гипотезу? Нет
N(0,1), 1000 наблюдений
statistic = 0.1886736348968725
Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений
statistic = 0.7096445710908483
Есть ли основания отбросить гипотезу? Нет
R:
%%R
res <- nortest::ad.test(data1)
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- nortest::ad.test(data2)
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- nortest::ad.test(data3)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- nortest::ad.test(data4)
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------ Anderson-Darling normality test data: data1 A = 0.29819, p-value = 0.5812 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,2), 100 наблюдений ------------------ Anderson-Darling normality test data: data2 A = 0.21992, p-value = 0.8309 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,1), 1000 наблюдений ------------------ Anderson-Darling normality test data: data3 A = 0.25351, p-value = 0.7325 Есть ли основания отбросить гипотезу? Нет ------------------ N(2,1), 1000 наблюдений ------------------ Anderson-Darling normality test data: data4 A = 0.59499, p-value = 0.1201 Есть ли основания отбросить гипотезу? Нет
т.е.данные получены из распределения с заданными параметрами.
Положим $F(x, \theta_0) \equiv \Phi(x)$, чтобы проверить, является ли распределение стандартным нормальным.
Python:
res = cramervonmises(data1, 'norm')
print('N(0,1), 100 наблюдений:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < 0.05) else 'Нет')
res = cramervonmises(data2, 'norm')
print('N(0,2), 100 наблюдений:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < 0.05) else 'Нет')
res = cramervonmises(data3, 'norm')
print('N(0,1), 1000 наблюдений:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < 0.05) else 'Нет')
res = cramervonmises(data4, 'norm')
print('N(2,1), 1000 наблюдений:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
pvalue = 0.36766184664388846
Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
pvalue = 0.0007999760848421689
Есть ли основания отбросить гипотезу? Да
N(0,1), 1000 наблюдений:
pvalue = 0.6078851470500846
Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
pvalue = 5.3411167155736905e-08
Есть ли основания отбросить гипотезу? Да
R:
%%R
res <- goftest::cvm.test(data1, "pnorm")
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- goftest::cvm.test(data2, "pnorm")
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- goftest::cvm.test(data3, "pnorm")
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- goftest::cvm.test(data4, "pnorm")
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution Parameters assumed to be fixed data: data1 omega2 = 0.080434, p-value = 0.6904 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,2), 100 наблюдений ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution Parameters assumed to be fixed data: data2 omega2 = 1.786, p-value = 3.433e-05 Есть ли основания отбросить гипотезу? Да ------------------ N(0,1), 1000 наблюдений ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution Parameters assumed to be fixed data: data3 omega2 = 0.032382, p-value = 0.9678 Есть ли основания отбросить гипотезу? Нет ------------------ N(2,1), 1000 наблюдений ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution Parameters assumed to be fixed data: data4 omega2 = 226.87, p-value < 2.2e-16 Есть ли основания отбросить гипотезу? Да
т.е. выборка имеет нормальное распределение, параметры которого заранее не известны.
Python:
_, pvalue = lilliefors(data1, 'norm')
print('N(0,1), 100 наблюдений:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
_, pvalue = lilliefors(data2, 'norm')
print('N(0,2), 100 наблюдений:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
_, pvalue = lilliefors(data3, 'norm')
print('N(0,1), 1000 наблюдений:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
_, pvalue = lilliefors(data4, 'norm')
print('N(2,1), 1000 наблюдений:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
pvalue = 0.5508216991827704
Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
pvalue = 0.28541751177174657
Есть ли основания отбросить гипотезу? Нет
N(0,1), 1000 наблюдений:
pvalue = 0.7354643954471619
Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
pvalue = 0.03457540982170773
Есть ли основания отбросить гипотезу? Да
R:
%%R
res <- lillie.test(data1)
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- lillie.test(data2)
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- lillie.test(data3)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- lillie.test(data4)
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: data1 D = 0.047378, p-value = 0.8401 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,2), 100 наблюдений ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: data2 D = 0.050275, p-value = 0.7727 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,1), 1000 наблюдений ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: data3 D = 0.017039, p-value = 0.6861 Есть ли основания отбросить гипотезу? Нет ------------------ N(2,1), 1000 наблюдений ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: data4 D = 0.026296, p-value = 0.09575 Есть ли основания отбросить гипотезу? Нет
т.е. выборка имеет нормальное распределение, параметры которого заранее не известны.
Критерий Шапиро - Франсия - это упрощение критерия Шапиро - Уилка и для него тоже необходимо, чтобы размер выборки был не более 5000.
Python:
res = shapiroFrancia(data1)
print('N(0,1), 100 наблюдений:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < 0.05) else 'Нет')
res = shapiroFrancia(data2)
print('N(0,2), 100 наблюдений:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < 0.05) else 'Нет')
res = shapiroFrancia(data3)
print('N(0,1), 1000 наблюдений:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < 0.05) else 'Нет')
res = shapiroFrancia(data4)
print('N(2,1), 1000 наблюдений:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
pvalue = 0.5539857967110066
Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
pvalue = 0.030175322690098563
Есть ли основания отбросить гипотезу? Да
N(0,1), 1000 наблюдений:
pvalue = 0.9786190291350362
Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
pvalue = 0.27727193178255205
Есть ли основания отбросить гипотезу? Нет
R:
%%R
res <- sf.test(data1)
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- sf.test(data2)
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- sf.test(data3)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
res <- sf.test(data4)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------ Shapiro-Francia normality test data: data1 W = 0.98897, p-value = 0.4974 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,2), 100 наблюдений ------------------ Shapiro-Francia normality test data: data2 W = 0.99344, p-value = 0.8462 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,1), 1000 наблюдений ------------------ Shapiro-Francia normality test data: data3 W = 0.99806, p-value = 0.2742 Есть ли основания отбросить гипотезу? Нет ------------------ N(0,1), 1000 наблюдений ------------------ Shapiro-Francia normality test data: data4 W = 0.99709, p-value = 0.06264 Есть ли основания отбросить гипотезу? Нет
Теперь исследуем теми же методами реальные данные. Рассмотрим следующие выборки:
Заметим, что выборки для драм и комедий пересекаются всего по 2 измерениям (см. код ниже).
Python:
# общие отрезки разбиения для гистограмм age
edges = np.histogram_bin_edges(df1['age'])
plt.figure(figsize=(10, 6))
sns.histplot(data=df1['age'], stat='density', bins=edges, kde=True)
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('age', fontsize=10, fontweight='bold')
plt.show()
print('Жен:', df1[df1.gender == 'F'].shape[0])
print('Муж:', df1[df1.gender == 'M'].shape[0])
Жен: 58 Муж: 117
plt.figure(figsize=(10, 6))
sns.histplot(data=df1[df1.gender == 'M'].age, stat='density', kde=True,
bins=edges, alpha=0.3, label='Male', color='blue')
sns.histplot(data=df1[df1.gender == 'F'].age, stat='density', kde=True,
bins=edges, alpha=0.3, label='Female', color='red')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('age (жен и муж)', fontsize=10, fontweight='bold')
plt.legend(title='Пол', title_fontsize=10, fontsize=10)
plt.show()
# общие отрезки разбиения для гистограмм average_rating
# (количество такое, так как average_rating округляется до 2 знаков после запятой)
edges = np.histogram_bin_edges(df2['average_rating'], bins=100)
plt.figure(figsize=(10, 6))
sns.histplot(data=df2['average_rating'], stat='density', bins=edges, kde=True)
plt.xticks(fontsize=9)
plt.xlabel('Средняя оценка', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('average_rating', fontsize=10, fontweight='bold')
plt.show()
df2_explode = df2.explode('genres')
comedy_df = df2_explode[df2_explode.genres == 'Comedy']
print('Комедии:', comedy_df.shape[0])
drama_df = df2_explode[df2_explode.genres == 'Drama']
print('Драмы:', drama_df.shape[0])
dramedy_df = pd.merge(comedy_df, drama_df, how='inner', on=['title'])
print('Драмеди:', dramedy_df.shape[0])
Комедии: 145 Драмы: 406 Драмеди: 2
plt.figure(figsize=(10, 6))
sns.histplot(data=comedy_df['average_rating'], stat='density', kde=True,
bins=edges, alpha=0.3, label='Comedy')
sns.histplot(data=drama_df['average_rating'], stat='density', kde=True,
bins=edges, alpha=0.3, label='Drama')
plt.xticks(fontsize=9)
plt.xlabel('Средняя оценка', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('average_rating (по жанрам)', fontsize=10, fontweight='bold')
plt.legend(title='Жанр', title_fontsize=10, fontsize=10)
plt.show()
R:
%%R
# общие отрезки разбиения для гистограмм age
edges = hist(df1$age, plot = FALSE)
%%R -w 1000 -h 600
hist(df1$age, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Возраст", ylab = "",
main = "age")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 2)
%%R
cat("Жен:", nrow(df1[df1$gender %like% "F", ]), "\n")
cat("Муж:", nrow(df1[df1$gender %like% "M", ]), "\n")
Жен: 58 Муж: 117
%%R -w 1000 -h 600
hist(df1[df1$gender %like% "F", "age"], freq = FALSE,
breaks = edges$breaks,
col = my_red,
xlab = "Возраст", ylab = "",
main = "age (муж и жен)")
lines(density(df1[df1$gender %like% "F", "age"], adjust=), col = "red", lwd = 3)
hist(df1[df1$gender %like% "M", "age"], freq = FALSE,
breaks = edges$breaks,
col = my_blue,
add = TRUE)
lines(density(df1[df1$gender %like% "M", "age"], adjust=), col = "royalblue", lwd = 3)
legend("topleft", inset=c(0.03, 0), legend=c("Male", "Female"), fill = c(my_blue, my_red), title="Пол")
%%R
# общие отрезки разбиения для гистограмм average_rating
# (количество такое, так как average_rating округляется до 2 знаков после запятой)
edges = hist(df2$average_rating, freq = FALSE, breaks = 100, plot = FALSE)
%%R -w 1000 -h 600
hist(df2$average_rating, freq = FALSE, breaks = edges$breaks,
col = "lightblue",
xlab = "Средняя оценка", ylab = "",
main = "average_rating")
lines(density(df2$average_rating, adjust=), col = "royalblue", lwd = 2)
%%R
comedy_df <- df2[df2$genres %like% "Comedy", ]
cat("Комедии:", nrow(comedy_df), "\n")
drama_df <- df2[df2$genres %like% "Drama", ]
cat("Драмы:", nrow(drama_df), "\n")
dramedy_df <- merge(comedy_df, drama_df, by = "title")
cat("Драмеди:", nrow(dramedy_df), "\n")
Комедии: 145 Драмы: 406 Драмеди: 2
%%R -w 1000 -h 600
my_orange <- t_col("orange", perc = 60)
hist(drama_df$average_rating, freq = FALSE,
breaks = edges$breaks,
col = my_orange,
xlab = "Средняя оценка", ylab = "",
main = "average_rating (по жанрам)")
lines(density(drama_df$average_rating, adjust=), col = "orange", lwd = 3)
hist(comedy_df$average_rating, freq = FALSE,
breaks = edges$breaks,
col = my_blue,
add = TRUE)
lines(density(comedy_df$average_rating, adjust=), col = "royalblue", lwd = 3)
legend("topleft", inset=c(0.03, 0), legend=c("Comedy", "Drama"), fill = c(my_blue, my_orange), title="Жанр")
Будем сравнивать данные с нормальными распределениями, которые имеют соответствующие средние и стандартные отклонения.
Python:
age_mean = statistics.mean(df1.age)
age_std = np.std(df1.age)
female_mean = statistics.mean(df1[df1.gender == 'F'].age)
female_std = np.std(df1[df1.gender == 'F'].age)
male_mean = statistics.mean(df1[df1.gender == 'M'].age)
male_std = np.std(df1[df1.gender == 'M'].age)
fig = plt.figure(figsize=(10,9))
x = np.linspace(0, 100, 1000)
ax1 = plt.subplot(221)
plt.plot(x, stats.norm.cdf(x, age_mean, age_std), color='red')
sns.ecdfplot(data=df1.age)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax1.set_title('age (все)', fontsize=10, fontweight='bold')
ax3 = plt.subplot(223)
plt.plot(x, stats.norm.cdf(x, male_mean, male_std), color='red')
sns.ecdfplot(data=df1[df1.gender == 'F'].age)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax3.set_title('age (жен)', fontsize=10, fontweight='bold')
ax4 = plt.subplot(224)
plt.plot(x, stats.norm.cdf(x, female_mean, female_std), color='red')
sns.ecdfplot(data=df1[df1.gender == 'M'].age)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax4.set_title('age (муж)', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=2)
plt.show()
all_mean = statistics.mean(df2.average_rating)
all_std = np.std(df2.average_rating)
comedy_mean = statistics.mean(comedy_df.average_rating)
comedy_std = np.std(comedy_df.average_rating)
drama_mean = statistics.mean(drama_df.average_rating)
drama_std = np.std(drama_df.average_rating)
fig = plt.figure(figsize=(10,9))
x = np.linspace(0, 8, 1000)
ax1 = plt.subplot(221)
plt.plot(x, stats.norm.cdf(x, all_mean, all_std), color='red')
sns.ecdfplot(data=df2.average_rating)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax1.set_title('average_rating (все)', fontsize=10, fontweight='bold')
ax3 = plt.subplot(223)
plt.plot(x, stats.norm.cdf(x, comedy_mean, comedy_std), color='red')
sns.ecdfplot(data=comedy_df.average_rating)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax3.set_title('average_rating (комедии)', fontsize=10, fontweight='bold')
ax4 = plt.subplot(224)
plt.plot(x, stats.norm.cdf(x, drama_mean, drama_std), color='red')
sns.ecdfplot(data=drama_df.average_rating)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax4.set_title('average_rating (драмы)', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=2)
plt.show()
R:
%%R
age_mean = mean(df1$age)
age_std = sd(df1$age)
female_mean = mean(df1[df1$gender %like% "F", "age"])
female_std = sd(df1[df1$gender %like% "F", "age"])
male_mean = mean(df1[df1$gender %like% "M", "age"])
male_std = sd(df1[df1$gender %like% "M", "age"])
%%R -w 1000 -h 450
plot(sort(df1$age), pnorm(sort(df1$age), mean = age_mean, sd = age_std), type = "l",
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="age (все)")
plot(ecdf(df1$age), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
par(mfrow = c(1,2))
plot(sort(df1[df1$gender %like% "F", "age"]), pnorm(sort(df1[df1$gender %like% "F", "age"]), mean = female_mean, sd = female_std), type = "l",
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="age (жен)")
plot(ecdf(df1[df1$gender %like% "F", "age"]), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
plot(sort(df1[df1$gender %like% "M", "age"]), pnorm(sort(df1[df1$gender %like% "M", "age"]), mean = male_mean, sd = male_std), type = "l",
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="age (муж)")
plot(ecdf(df1[df1$gender %like% "M", "age"]), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
%%R
all_mean = mean(df2$average_rating)
all_std = sd(df2$average_rating)
comedy_mean = mean(comedy_df$average_rating)
comedy_std = sd(comedy_df$average_rating)
drama_mean = mean(drama_df$average_rating)
drama_std = sd(drama_df$average_rating)
%%R -w 1000 -h 450
plot(sort(df2$average_rating), pnorm(sort(df2$average_rating), mean = all_mean, sd = all_std), type = "l",
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="average_rating (все)")
plot(ecdf(df2$average_rating), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
par(mfrow = c(1,2))
plot(sort(comedy_df$average_rating), pnorm(sort(comedy_df$average_rating), mean = comedy_mean, sd = comedy_std), type = "l",
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="average_rating (комедии)")
plot(ecdf(comedy_df$average_rating), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
plot(sort(drama_df$average_rating), pnorm(sort(drama_df$average_rating), mean = drama_mean, sd = drama_std), type = "l",
lwd = 2, col = "red",
xlab = "x", ylab = "F(x)",
main="average_rating (драмы)")
plot(ecdf(drama_df$average_rating), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
Python:
fig = plt.figure(figsize=(10,9))
ax1 = plt.subplot(221)
sm.qqplot(df1.age, ax=ax1, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax1.set_title('age (все)', fontsize=10, fontweight='bold')
ax3 = plt.subplot(223)
sm.qqplot(df1[df1.gender == 'F'].age, ax=ax3, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax3.set_title('age (жен)', fontsize=10, fontweight='bold')
ax4 = plt.subplot(224)
sm.qqplot(df1[df1.gender == 'M'].age, ax=ax4, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax4.set_title('age (муж)', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=2)
plt.show()
fig = plt.figure(figsize=(10,9))
ax1 = plt.subplot(221)
sm.qqplot(df2.average_rating, ax=ax1, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax1.set_title('average_rating (все)', fontsize=10, fontweight='bold')
ax3 = plt.subplot(223)
sm.qqplot(comedy_df.average_rating, ax=ax3, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax3.set_title('average_rating (комедии)', fontsize=10, fontweight='bold')
ax4 = plt.subplot(224)
sm.qqplot(drama_df.average_rating, ax=ax4, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax4.set_title('average_rating (драмы)', fontsize=10, fontweight='bold')
fig.tight_layout(h_pad=2, w_pad=2)
plt.show()
R:
%%R -w 1000 -h 450
qqnorm(df1$age,
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "age (все)")
qqline(df1$age,
col = "red", lwd = 2)
par(mfrow = c(1,2))
qqnorm(df1[df1$gender %like% "F", "age"],
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "age (жен)")
qqline(df1[df1$gender %like% "F", "age"],
col = "red", lwd = 2)
qqnorm(df1[df1$gender %like% "M", "age"],
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "age (муж)")
qqline(df1[df1$gender %like% "M", "age"],
col = "red", lwd = 2)
%%R -w 1000 -h 450
qqnorm(df2$average_rating,
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "average_rating (все)")
qqline(df2$average_rating,
col = "red", lwd = 2)
par(mfrow = c(1,2))
qqnorm(comedy_df$average_rating,
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "average_rating (комедии)")
qqline(comedy_df$average_rating,
col = "red", lwd = 2)
qqnorm(drama_df$average_rating,
col = "royalblue", pch = 19,
xlab = "Теоретические квантили", ylab = "Выборочные квантили",
main = "average_rating (драмы)")
qqline(drama_df$average_rating,
col = "red", lwd = 2)
Python:
fig, ax = plt.subplots(nrows=2, ncols=2, figsize = (10, 9))
envelmet(df1.age, 'age (все)', ax[0, 0])
fig.delaxes(ax[0, 1])
envelmet(df1[df1.gender == 'F'].age, 'age (жен)', ax[1, 0])
envelmet(df1[df1.gender == 'M'].age, 'age (муж)', ax[1, 1])
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(nrows=2, ncols=2, figsize = (10, 9))
envelmet(df2.average_rating, 'average_rating (все)', ax[0, 0])
fig.delaxes(ax[0, 1])
envelmet(comedy_df.average_rating, 'average_rating (комедии)', ax[1, 0])
envelmet(drama_df.average_rating, 'average_rating (драмы)', ax[1, 1])
fig.tight_layout()
plt.show()
R:
%%R -w 1000 -h 450
envelmet(df1$age, "age (все)")
par(mfrow = c(1,2))
envelmet(df1[df1$gender %like% "F", "age"], "age (жен)")
envelmet(df1[df1$gender %like% "M", "age"], "age (муж)")
%%R -w 1000 -h 450
envelmet(df2$average_rating, "average_rating (все)")
par(mfrow = c(1,2))
envelmet(comedy_df$average_rating, "average_rating (комедии)")
envelmet(drama_df$average_rating, "average_rating (драмы)")
Воспользуемся двухвыборочным критерием.
Есть 2 выборки $X_1, \dots, X_n$ и $Y_1, \dots, Y_m$ с функциями распределения $F(x)$ и $G(x)$ соответственно. $$H_0: F(x) \equiv G(x),$$ т.е. выборки принадлежат к одному распределению.
Если полученный при использовании критерия $pvalue$ меньше заданного заранее уровня значимости $\alpha$ (положим $\alpha = 0.5$), то $H_0$ должна быть отвергнута, иначе - нет оснований отбросить нулевую гипотезу.
Во-первых, чтобы проверить нормальность данных, семплируем по выборке из нормальных распределений со средними и стандартными отклонениями, соответствующими этим же величинам у исследуемых распределений.
Во-вторых, проверим, одинаково ли распределены данные:
Python:
alpha = 0.05
age_norm = np.random.normal(age_mean, age_std, size=1000)
female_norm = np.random.normal(female_mean, female_std, size=100)
male_norm = np.random.normal(male_mean, male_mean, size=100)
_, pvalue = kstest(df1.age, age_norm)
print('age - нормальное распределение:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
_, pvalue = kstest(df1[df1.gender == 'F'].age, female_norm)
print('age (жен) - нормальное распределение:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
_, pvalue = kstest(df1[df1.gender == 'M'].age, male_norm)
print('age (муж) - нормальное распределение:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
age - нормальное распределение:
pvalue = 0.12182158757492134
Есть ли основания отбросить гипотезу? Нет
age (жен) - нормальное распределение:
pvalue = 0.910886686614063
Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
pvalue = 6.907798016052148e-12
Есть ли основания отбросить гипотезу? Да
_, pvalue = kstest(df1[df1.gender == 'F'].age, df1[df1.gender == 'M'].age)
print('Данные по женщинам и мужчинам одинаково распределены:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
Данные по женщинам и мужчинам одинаково распределены:
pvalue = 0.7227115888214679
Есть ли основания отбросить гипотезу? Нет
all_norm = np.random.normal(all_mean, all_std, size=1000)
comedy_norm = np.random.normal(comedy_mean, comedy_std, size=100)
drama_norm = np.random.normal(drama_mean, comedy_mean, size=100)
_, pvalue = kstest(df2.average_rating, all_norm)
print('average_rating - нормальное распределение:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
_, pvalue = kstest(comedy_df.average_rating, comedy_norm)
print('average_rating (комедии) - нормальное распределение:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
_, pvalue = kstest(drama_df.average_rating, drama_norm)
print('average_rating (драмы) - нормальное распределение:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
average_rating - нормальное распределение:
pvalue = 0.001585106579657839
Есть ли основания отбросить гипотезу? Да
average_rating (комедии) - нормальное распределение:
pvalue = 0.8805561402063771
Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
pvalue = 8.713468542532803e-14
Есть ли основания отбросить гипотезу? Да
_, pvalue = kstest(comedy_df.average_rating, drama_df.average_rating)
print('Данные по комедиям и драмам одинаково распределены:')
print(' pvalue =', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
Данные по комедиям и драмам одинаково распределены:
pvalue = 0.020898085404186383
Есть ли основания отбросить гипотезу? Да
R:
%%R
alpha <- 0.05
%%R
age_norm <- rnorm(1000, mean = age_mean, sd = age_std)
female_norm <- rnorm(100, mean = female_mean, sd = female_std)
male_norm <- rnorm(100, mean = male_mean, sd = male_mean)
%%R
res <- ks.test(df1$age, age_norm)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- ks.test(df1[df1$gender %like% "F", "age"], female_norm)
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- ks.test(df1[df1$gender %like% "M", "age"], male_norm)
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------ Asymptotic two-sample Kolmogorov-Smirnov test data: df1$age and age_norm D = 0.11271, p-value = 0.04545 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Да ------------------ age (жен) - нормальное распределение ------------------ Exact two-sample Kolmogorov-Smirnov test data: df1[df1$gender %like% "F", "age"] and female_norm D = 0.20379, p-value = 0.07555 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Нет ------------------ age (муж) - нормальное распределение ------------------ Asymptotic two-sample Kolmogorov-Smirnov test data: df1[df1$gender %like% "M", "age"] and male_norm D = 0.42291, p-value = 8.418e-09 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Да
%%R
res <- ks.test(df1[df1$gender %like% "F", "age"], df1[df1$gender %like% "M", "age"])
cat("------------------ Данные по женщинам и мужчинам одинаково распределены ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ Данные по женщинам и мужчинам одинаково распределены ------------------ Exact two-sample Kolmogorov-Smirnov test data: df1[df1$gender %like% "F", "age"] and df1[df1$gender %like% "M", "age"] D = 0.1064, p-value = 0.6048 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Нет
%%R
all_norm <- rnorm(1000, mean = all_mean, sd = all_std)
comedy_norm <- rnorm(100, mean = comedy_mean, sd = comedy_std)
drama_norm <- rnorm(100, mean = drama_mean, sd = comedy_mean)
%%R
res <- ks.test(df2$average_rating, all_norm)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- ks.test(comedy_df$average_rating, comedy_norm)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- ks.test(drama_df$average_rating, drama_norm)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------ Asymptotic two-sample Kolmogorov-Smirnov test data: df2$average_rating and all_norm D = 0.0667, p-value = 0.0006139 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Да ------------------ average_rating (комедии) - нормальное распределение ------------------ Asymptotic two-sample Kolmogorov-Smirnov test data: comedy_df$average_rating and comedy_norm D = 0.10724, p-value = 0.504 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Нет ------------------ average_rating (драмы) - нормальное распределение ------------------ Asymptotic two-sample Kolmogorov-Smirnov test data: drama_df$average_rating and drama_norm D = 0.43015, p-value = 2.547e-13 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Да
%%R
res <- ks.test(comedy_df$average_rating, drama_df$average_rating)
cat("------------------ Данные по комедиям и драмам одинаково распределены ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ Данные по комедиям и драмам одинаково распределены ------------------ Asymptotic two-sample Kolmogorov-Smirnov test data: comedy_df$average_rating and drama_df$average_rating D = 0.14433, p-value = 0.02332 alternative hypothesis: two-sided Есть ли основания отбросить гипотезу? Да
Так как для критерия Шапиро-Уилка необходимо, чтобы размер выборки был не более 5000, возьмём только половину от всех данных по рейтингам (отсортируем по возрастанию average_rating и возьмём каждое второе значение).
Python:
_, pvalue = shapiro(df1.age)
print('age - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, pvalue = shapiro(df1[df1.gender == 'F'].age)
print('age (жен) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, pvalue = shapiro(df1[df1.gender == 'M'].age)
print('age (муж) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
age - нормальное распределение:
pvalue = 0.0003017295675817877
Есть ли основания отбросить гипотезу? Да
age (жен) - нормальное распределение:
pvalue = 0.24597471952438354
Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
pvalue = 0.0003075910790357739
Есть ли основания отбросить гипотезу? Да
df2_half = df2.sort_values('average_rating').reset_index().iloc[lambda x: x.index % 2 == 0].reset_index()
_, pvalue = shapiro(df2_half.average_rating)
print('average_rating - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, pvalue = shapiro(comedy_df.average_rating)
print('average_rating (комедии) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, pvalue = shapiro(drama_df.average_rating)
print('average_rating (драмы) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
average_rating - нормальное распределение:
pvalue = 0.0
Есть ли основания отбросить гипотезу? Да
average_rating (комедии) - нормальное распределение:
pvalue = 0.7259739637374878
Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
pvalue = 3.4041782726035308e-09
Есть ли основания отбросить гипотезу? Да
R:
%%R
res <- shapiro.test(df1$age)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- shapiro.test(df1[df1$gender %like% "F", "age"])
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- shapiro.test(df1[df1$gender %like% "M", "age"])
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------ Shapiro-Wilk normality test data: df1$age W = 0.96624, p-value = 0.0003017 Есть ли основания отбросить гипотезу? Да ------------------ age (жен) - нормальное распределение ------------------ Shapiro-Wilk normality test data: df1[df1$gender %like% "F", "age"] W = 0.97398, p-value = 0.246 Есть ли основания отбросить гипотезу? Нет ------------------ age (муж) - нормальное распределение ------------------ Shapiro-Wilk normality test data: df1[df1$gender %like% "M", "age"] W = 0.95093, p-value = 0.0003076 Есть ли основания отбросить гипотезу? Да
%%R
df2_half <- df2 %>%
arrange(average_rating) %>%
mutate(index = row_number()) %>%
filter(index %% 2 == 0) %>%
select(-index) %>%
mutate(index = row_number()) %>%
select(-index)
%%R
res <- shapiro.test(df2_half$average_rating)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- shapiro.test(comedy_df$average_rating)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- shapiro.test(drama_df$average_rating)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------ Shapiro-Wilk normality test data: df2_half$average_rating W = 0.91881, p-value < 2.2e-16 Есть ли основания отбросить гипотезу? Да ------------------ average_rating (комедии) - нормальное распределение ------------------ Shapiro-Wilk normality test data: comedy_df$average_rating W = 0.99322, p-value = 0.726 Есть ли основания отбросить гипотезу? Нет ------------------ average_rating (драмы) - нормальное распределение ------------------ Shapiro-Wilk normality test data: drama_df$average_rating W = 0.95912, p-value = 3.405e-09 Есть ли основания отбросить гипотезу? Да
Python:
На Python есть реализация критерия Андерсона для $k$ выборок. С его помощью (также как с двухвыборочным критерием Колмогорова - Смирнова) проверим и то, что данные нормальны, и то, что данные по комедиям и драмам распределены одинаковы.
_, _, pvalue = anderson_ksamp([df1.age, age_norm])
print('age - нормальное распределеныие:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, _, pvalue = anderson_ksamp([df1[df1.gender == 'F'].age, female_norm])
print('age (жен) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, _, pvalue = anderson_ksamp([df1[df1.gender == 'M'].age, male_norm])
print('age (муж) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
age - нормальное распределеныие:
pvalue = 0.25
Есть ли основания отбросить гипотезу? Нет
age (жен) - нормальное распределение:
pvalue = 0.25
Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
pvalue = 0.001
Есть ли основания отбросить гипотезу? Да
_, _, pvalue = anderson_ksamp([df1[df1.gender == 'F'].age, df1[df1.gender == 'M'].age])
print('Данные по женщинам и мужчинам одинаково распределены:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
Данные по женщинам и мужчинам одинаково распределены:
pvalue = 0.25
Есть ли основания отбросить гипотезу? Нет
_, _, pvalue = anderson_ksamp([df2.average_rating, all_norm])
print('average_rating - нормальное распределеныие:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, _, pvalue = anderson_ksamp([comedy_df.average_rating, comedy_norm])
print('average_rating (комедии) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, _, pvalue = anderson_ksamp([drama_df.average_rating, drama_norm])
print('average_rating (драмы) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
average_rating - нормальное распределеныие:
pvalue = 0.0033431320434570704
Есть ли основания отбросить гипотезу? Да
average_rating (комедии) - нормальное распределение:
pvalue = 0.25
Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
pvalue = 0.001
Есть ли основания отбросить гипотезу? Да
_, _, pvalue = anderson_ksamp([comedy_df.average_rating, drama_df.average_rating])
print('Данные по комедиям и драмам одинаково распределены:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
Данные по комедиям и драмам одинаково распределены:
pvalue = 0.03658924476599451
Есть ли основания отбросить гипотезу? Да
R:
%%R
res <- nortest::ad.test(df1$age)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- nortest::ad.test(df1[df1$gender %like% "F", "age"])
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- nortest::ad.test(df1[df1$gender %like% "M", "age"])
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------ Anderson-Darling normality test data: df1$age A = 1.3089, p-value = 0.002072 Есть ли основания отбросить гипотезу? Да ------------------ age (жен) - нормальное распределение ------------------ Anderson-Darling normality test data: df1[df1$gender %like% "F", "age"] A = 0.5091, p-value = 0.1906 Есть ли основания отбросить гипотезу? Нет ------------------ age (муж) - нормальное распределение ------------------ Anderson-Darling normality test data: df1[df1$gender %like% "M", "age"] A = 1.0999, p-value = 0.006709 Есть ли основания отбросить гипотезу? Да
%%R
res <- nortest::ad.test(df2$average_rating)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- nortest::ad.test(comedy_df$average_rating)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- nortest::ad.test(drama_df$average_rating)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------ Anderson-Darling normality test data: df2$average_rating A = 57.075, p-value < 2.2e-16 Есть ли основания отбросить гипотезу? Да ------------------ average_rating (комедии) - нормальное распределение ------------------ Anderson-Darling normality test data: comedy_df$average_rating A = 0.30901, p-value = 0.5541 Есть ли основания отбросить гипотезу? Нет ------------------ average_rating (драмы) - нормальное распределение ------------------ Anderson-Darling normality test data: drama_df$average_rating A = 4.5655, p-value = 2.455e-11 Есть ли основания отбросить гипотезу? Да
Python:
На Python есть реализация критерия Крамера - фон Мизеса для $k$ выборок. С его помощью (также как и с двухвыборочным критерием Колмогорова - Смирнова) проверим и то, что данные нормальны, и то, что данные по мужчинам/женщинам и комедиям/драмам распределены одинаковы.
res = cramervonmises_2samp(df1.age, age_norm)
print('age - нормальное распределение:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
res = cramervonmises_2samp(df1[df1.gender == 'F'].age, female_norm)
print('age (жен) - нормальное распределение:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
res = cramervonmises_2samp(df1[df1.gender == 'M'].age, male_norm)
print('age (муж) - нормальное распределение:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
age - нормальное распределение:
pvalue = 0.39133160782598164
Есть ли основания отбросить гипотезу? Нет
age (жен) - нормальное распределение:
pvalue = 0.9345779001744031
Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
pvalue = 1.1599448268562185e-07
Есть ли основания отбросить гипотезу? Да
res = cramervonmises_2samp(df1[df1.gender == 'F'].age, df1[df1.gender == 'M'].age)
print('Данные по женщинам и мужчинам одинаково распределены:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
Данные по женщинам и мужчинам одинаково распределены:
pvalue = 0.5198899109819444
Есть ли основания отбросить гипотезу? Нет
res = cramervonmises_2samp(df2.average_rating, all_norm)
print('average_rating - нормальное распределение:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
res = cramervonmises_2samp(comedy_df.average_rating, comedy_norm)
print('average_rating (комедии) - нормальное распределение:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
res = cramervonmises_2samp(drama_df.average_rating, drama_norm)
print('average_rating (драмы) - нормальное распределение:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
average_rating - нормальное распределение:
pvalue = 0.00012693698164067957
Есть ли основания отбросить гипотезу? Да
average_rating (комедии) - нормальное распределение:
pvalue = 0.8909349351026549
Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
pvalue = 2.1094737068239056e-10
Есть ли основания отбросить гипотезу? Да
res = cramervonmises_2samp(comedy_df.average_rating, drama_df.average_rating)
print('Данные по комедиям и драмам одинаково распределены:')
print(' pvalue = ', res.pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
Данные по комедиям и драмам одинаково распределены:
pvalue = 0.03926253837913829
Есть ли основания отбросить гипотезу? Да
R:
%%R
res <- goftest::cvm.test(df1$age, "pnorm", mean = age_mean, sd = age_std)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- goftest::cvm.test(df1[df1$gender %like% 'F', "age"], "pnorm", mean = female_mean, sd = female_std)
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- goftest::cvm.test(df1[df1$gender %like% 'M', "age"], "pnorm", mean = male_mean, sd = male_std)
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution with parameters mean = 42.8457142857143, sd = 10.9052434624528 Parameters assumed to be fixed data: df1$age omega2 = 0.19279, p-value = 0.282 Есть ли основания отбросить гипотезу? Нет ------------------ age (жен) - нормальное распределение ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution with parameters mean = 41.6206896551724, sd = 10.8786408887419 Parameters assumed to be fixed data: df1[df1$gender %like% "F", "age"] omega2 = 0.084203, p-value = 0.6697 Есть ли основания отбросить гипотезу? Нет ------------------ age (муж) - нормальное распределение ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution with parameters mean = 43.4529914529915, sd = 10.913844202865 Parameters assumed to be fixed data: df1[df1$gender %like% "M", "age"] omega2 = 0.14668, p-value = 0.4003 Есть ли основания отбросить гипотезу? Нет
%%R
res <- goftest::cvm.test(df2$average_rating, "pnorm", mean = all_mean, sd = all_std)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- goftest::cvm.test(comedy_df$average_rating, "pnorm", mean = comedy_mean, sd = comedy_std)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- goftest::cvm.test(drama_df$average_rating, "pnorm", mean = drama_mean, sd = drama_std)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution with parameters mean = 4.068577, sd = 0.335358601515096 Parameters assumed to be fixed data: df2$average_rating omega2 = 8.5192, p-value < 2.2e-16 Есть ли основания отбросить гипотезу? Да ------------------ average_rating (комедии) - нормальное распределение ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution with parameters mean = 4.02427586206897, sd = 0.266092467291777 Parameters assumed to be fixed data: comedy_df$average_rating omega2 = 0.038392, p-value = 0.9418 Есть ли основания отбросить гипотезу? Нет ------------------ average_rating (драмы) - нормальное распределение ------------------ Cramer-von Mises test of goodness-of-fit Null hypothesis: Normal distribution with parameters mean = 4.02997536945813, sd = 0.294512360107409 Parameters assumed to be fixed data: drama_df$average_rating omega2 = 0.80123, p-value = 0.007232 Есть ли основания отбросить гипотезу? Да
Python:
_, pvalue = lilliefors(df1.age, 'norm')
print('age - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, pvalue = lilliefors(df1[df1.gender == 'F'].age, 'norm')
print('age (жен) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, pvalue = lilliefors(df1[df1.gender == 'M'].age, 'norm')
print('age (муж) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
age - нормальное распределение:
pvalue = 0.001761024787129434
Есть ли основания отбросить гипотезу? Да
age (жен) - нормальное распределение:
pvalue = 0.26096030891637234
Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
pvalue = 0.0206531490309459
Есть ли основания отбросить гипотезу? Да
_, pvalue = lilliefors(df2.average_rating, 'norm')
print('average_rating - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, pvalue = lilliefors(comedy_df.average_rating, 'norm')
print('average_rating (комедии) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
_, pvalue = lilliefors(drama_df.average_rating, 'norm')
print('average_rating (драмы) - нормальное распределение:')
print(' pvalue = ', pvalue)
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
average_rating - нормальное распределение:
pvalue = 0.0009999999999998899
Есть ли основания отбросить гипотезу? Да
average_rating (комедии) - нормальное распределение:
pvalue = 0.9073222722620489
Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
pvalue = 0.0009999999999998899
Есть ли основания отбросить гипотезу? Да
R:
%%R
res <- lillie.test(df1$age)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- lillie.test(df1[df1$gender %like% 'F', "age"])
cat("------------------ age (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- lillie.test(df1[df1$gender %like% 'M', "age"])
cat("------------------ age (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: df1$age D = 0.092914, p-value = 0.0008559 Есть ли основания отбросить гипотезу? Да ------------------ age (комедии) - нормальное распределение ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: df1[df1$gender %like% "F", "age"] D = 0.093698, p-value = 0.2327 Есть ли основания отбросить гипотезу? Нет ------------------ age (драмы) - нормальное распределение ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: df1[df1$gender %like% "M", "age"] D = 0.09319, p-value = 0.01424 Есть ли основания отбросить гипотезу? Да
%%R
res <- lillie.test(df2$average_rating)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- lillie.test(comedy_df$average_rating)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- lillie.test(drama_df$average_rating)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: df2$average_rating D = 0.053433, p-value < 2.2e-16 Есть ли основания отбросить гипотезу? Да ------------------ average_rating (комедии) - нормальное распределение ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: comedy_df$average_rating D = 0.038619, p-value = 0.8609 Есть ли основания отбросить гипотезу? Нет ------------------ average_rating (драмы) - нормальное распределение ------------------ Lilliefors (Kolmogorov-Smirnov) normality test data: drama_df$average_rating D = 0.11116, p-value = 4.371e-13 Есть ли основания отбросить гипотезу? Да
Критерий Шапиро - Франсия - это упрощение критерия Шапиро - Уилка и для него тоже необходимо, чтобы размер выборки был не более 5000. Опять берём только половину данных по рейтингам.
Python:
res = shapiroFrancia(df1.first_published)
print('age - нормальное распределение:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')
res = shapiroFrancia(df1[df1.gender == 'F'].age)
print('age (жен) - нормальное распределение:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')
res = shapiroFrancia(df1[df1.gender == 'M'].age)
print('age (муж) - нормальное распределение:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')
age - нормальное распределение:
pvalue = 1.750186564963448e-18
Есть ли основания отбросить гипотезу? Да
age (жен) - нормальное распределение:
pvalue = 0.1360614634580179
Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
pvalue = 0.0005470998390365658
Есть ли основания отбросить гипотезу? Да
res = shapiroFrancia(df2_half.average_rating)
print('average_rating - нормальное распределение:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')
res = shapiroFrancia(comedy_df.average_rating)
print('average_rating (комедии) - нормальное распределение:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')
res = shapiroFrancia(drama_df.average_rating)
print('average_rating (драмы) - нормальное распределение:')
print(' pvalue = ', res['p-value'])
print(' Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')
average_rating - нормальное распределение:
pvalue = 3.293768000557143e-46
Есть ли основания отбросить гипотезу? Да
average_rating (комедии) - нормальное распределение:
pvalue = 0.46524972040289425
Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
pvalue = 2.0959308700711756e-08
Есть ли основания отбросить гипотезу? Да
R:
%%R
res <- sf.test(df1$age)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- sf.test(df1[df1$gender %like% "F", "age"])
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- sf.test(df1[df1$gender %like% "M", "age"])
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------ Shapiro-Francia normality test data: df1$age W = 0.96513, p-value = 0.0004398 Есть ли основания отбросить гипотезу? Да ------------------ age (жен) - нормальное распределение ------------------ Shapiro-Francia normality test data: df1[df1$gender %like% "F", "age"] W = 0.96964, p-value = 0.1361 Есть ли основания отбросить гипотезу? Нет ------------------ age (муж) - нормальное распределение ------------------ Shapiro-Francia normality test data: df1[df1$gender %like% "M", "age"] W = 0.95049, p-value = 0.0005471 Есть ли основания отбросить гипотезу? Да
%%R
res <- sf.test(df2_half$average_rating)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- sf.test(comedy_df$average_rating)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
res <- sf.test(drama_df$average_rating)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------ Shapiro-Francia normality test data: df2_half$average_rating W = 0.91803, p-value < 2.2e-16 Есть ли основания отбросить гипотезу? Да ------------------ average_rating (комедии) - нормальное распределение ------------------ Shapiro-Francia normality test data: comedy_df$average_rating W = 0.99158, p-value = 0.4652 Есть ли основания отбросить гипотезу? Нет ------------------ average_rating (драмы) - нормальное распределение ------------------ Shapiro-Francia normality test data: drama_df$average_rating W = 0.95811, p-value = 2.096e-08 Есть ли основания отбросить гипотезу? Да
Заметим:
В первом датасете:
При доверительном уровне $\alpha = 0.1$ картина не меняется. А при доверительном уровне $\alpha = 0.01$ данные по мужчинам проходят ещё и критерий Лиллиефорса.
Поэтому будем считать, что распределения данных по всем возрастам и по возрастам мужчин достаточно близки к нормальным, чтобы можно было пренебречь этой разницей (скорее всего, у этих распределений более тяжёлые хвосты: это видно из гистограммы; к тому же известно, что сильнее критерия Колмогорова - Смирнова критерий Крамера - фон Мизеса реагирует на отклонения в среднем, а критерий Андерсона - Дарлинга - на тяжёлые хвосты).
Во втором датасете:
Для каждого из критериев формулируется нулевая гипотеза $H_0$ и альтернативная гипотеза. Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается в пользу альтернативы, иначе - нет оснований отбросить гипотезу.
Также все разновидности критерия Стьюдента основаны на дополнительном предположении о нормальности выборки данных, поэтому будем использовать эти разновидности на данных о возрастах авторов, которые считаем нормально распределёнными.
Дана выборка $X_1, \dots, X_n \sim N(\mu, \alpha).$ $$H_0: \mathbb{E}X = \mu_0,$$ т.е. матожидание выборки равно заданному числу.
Альтернативные гипотезы: $$H_1: \mathbb{E}X \neq \mu_0 \quad (\text{двухсторонняя гипотеза})$$ $$H_1': \mathbb{E}X > \mu_0 \quad (\text{односторонняя гипотеза})$$ $$H_1'': \mathbb{E}X < \mu_0 \quad (\text{односторонняя гипотеза})$$
Проверим, равен ли средний возраст авторов 40 годам.
Python:
print('mean(X) =', statistics.mean(df1.age), 'vs mu0 = 40')
print('E(X) != mu0')
res = ttest(df1.age, 40, alternative = 'two-sided', confidence=0.9)
print(' 0.90: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'two-sided', confidence=0.95)
print(' 0.95: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'two-sided', confidence=0.99)
print(' 0.99: pvalue =', res['p-val'][0], ' power =', res['power'][0])
print('E(X) > mu0')
res = ttest(df1.age, 40, alternative = 'greater', confidence=0.9)
print(' 0.90: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'greater', confidence=0.95)
print(' 0.95: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'greater', confidence=0.99)
print(' 0.99: pvalue =', res['p-val'][0], ' power =', res['power'][0])
print('E(X) < mu0')
res = ttest(df1.age, 40, alternative = 'less', confidence=0.9)
print(' 0.90: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'less', confidence=0.95)
print(' 0.95: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'less', confidence=0.99)
print(' 0.99: pvalue =', res['p-val'][0], ' power =', res['power'][0])
mean(X) = 42.84571428571429 vs mu0 = 40
E(X) != mu0
0.90: pvalue = 0.0006984719380686418 power = 0.9296182976236684
0.95: pvalue = 0.0006984719380686418 power = 0.9296182976236684
0.99: pvalue = 0.0006984719380686418 power = 0.9296182976236684
E(X) > mu0
0.90: pvalue = 0.0003492359690343209 power = 0.9635704269241538
0.95: pvalue = 0.0003492359690343209 power = 0.9635704269241538
0.99: pvalue = 0.0003492359690343209 power = 0.9635704269241538
E(X) < mu0
0.90: pvalue = 0.9996507640309656 power = 1.8520793543252978e-07
0.95: pvalue = 0.9996507640309656 power = 1.8520793543252978e-07
0.99: pvalue = 0.9996507640309656 power = 1.8520793543252978e-07
R:
%%R
cat("mean(X) =", mean(df1$age), "vs mu0 = 40")
cat("\nE(X) == mu0\n")
res = t.test(df1$age, mu = 40, alternative = "two.sided", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "two.sided", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "two.sided", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, "\n")
cat("\nE(X) < mu0\n")
res = t.test(df1$age, mu = 40, alternative = "greater", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "greater", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "greater", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, "\n")
cat("\nE(X) > mu0\n")
res = t.test(df1$age, mu = 40, alternative = "less", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "less", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "less", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, "\n")
mean(X) = 42.84571 vs mu0 = 40
E(X) == mu0
0.90: pvalue = 0.0006984719
0.95: pvalue = 0.0006984719
0.99: pvalue = 0.0006984719
E(X) < mu0
0.90: pvalue = 0.000349236
0.95: pvalue = 0.000349236
0.99: pvalue = 0.000349236
E(X) > mu0
0.90: pvalue = 0.9996508
0.95: pvalue = 0.9996508
0.99: pvalue = 0.9996508
Средний возраст авторов больше 40.
Даны независимые выборки $X_1, \dots, X_n \sim N(\mu_1, \alpha_1)$ и $Y_1, \dots, Y_m \sim N(\mu_2, \alpha_2).$ $$H_0: \mathbb{E}X = \mathbb{E}Y,$$ т.е. матожидания выборок равны.
Альтернативные гипотезы: $$H_1: \mathbb{E}X \neq \mathbb{E}Y \quad (\text{двухсторонняя гипотеза})$$ $$H_1': \mathbb{E}X > \mathbb{E}Y \quad (\text{односторонняя гипотеза})$$ $$H_1'': \mathbb{E}X < \mathbb{E}Y \quad (\text{односторонняя гипотеза})$$
Проверим, равны ли средний возраст авторов-мужчин и авторов-женщин.
Python:
print('mean(X) =', statistics.mean(df1[df1.gender == 'M'].age),
'vs mean(Y) =', statistics.mean(df1[df1.gender == 'F'].age))
print('E(X) != E(Y)')
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'two-sided', confidence=0.9)
print(' 0.90: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'two-sided', confidence=0.95)
print(' 0.95: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'two-sided', confidence=0.99)
print(' 0.99: pvalue =', res['p-val'][0], ' power =', res['power'][0])
print('E(X) > E(Y)')
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'greater', confidence=0.9)
print(' 0.90: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'greater', confidence=0.95)
print(' 0.95: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'greater', confidence=0.99)
print(' 0.99: pvalue =', res['p-val'][0], ' power =', res['power'][0])
print('E(X) < E(Y)')
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'less', confidence=0.9)
print(' 0.90: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'less', confidence=0.95)
print(' 0.95: pvalue =', res['p-val'][0], ' power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'less', confidence=0.99)
print(' 0.99: pvalue =', res['p-val'][0], ' power =', res['power'][0])
mean(X) = 43.452991452991455 vs mean(Y) = 41.62068965517241
E(X) != E(Y)
0.90: pvalue = 0.2969828479314967 power = 0.1803430386782663
0.95: pvalue = 0.2969828479314967 power = 0.1803430386782663
0.99: pvalue = 0.2969828479314967 power = 0.1803430386782663
E(X) > E(Y)
0.90: pvalue = 0.14849142396574835 power = 0.2734624976802309
0.95: pvalue = 0.14849142396574835 power = 0.2734624976802309
0.99: pvalue = 0.14849142396574835 power = 0.2734624976802309
E(X) < E(Y)
0.90: pvalue = 0.8515085760342517 power = 0.0036011101618038895
0.95: pvalue = 0.8515085760342517 power = 0.0036011101618038895
0.99: pvalue = 0.8515085760342517 power = 0.0036011101618038895
R:
%%R
cat("mean(X) =", mean(df1[df1$gender %like% 'M', "age"]),
"vs mean(Y) =", mean(mean(df1[df1$gender %like% 'F', "age"])))
cat("\nE(X) != E(Y)\n")
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "two.sided", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "two.sided", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "two.sided", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
cat("\nE(X) > E(Y)\n")
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "greater", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "greater", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "greater", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
cat("\nE(X) < E(Y)\n")
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "less", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "less", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "less", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
mean(X) = 43.45299 vs mean(Y) = 41.62069
E(X) != E(Y)
0.90: pvalue = 0.2969828
0.95: pvalue = 0.2969828
0.99: pvalue = 0.2969828
E(X) > E(Y)
0.90: pvalue = 0.1484914
0.95: pvalue = 0.1484914
0.99: pvalue = 0.1484914
E(X) < E(Y)
0.90: pvalue = 0.8515086
0.95: pvalue = 0.8515086
0.99: pvalue = 0.8515086
Средний возраст авторов-мужчин выше среднего возраста авторов-женщин.
Даны независимые выборки $X_1, \dots, X_n$ и $Y_1, \dots, Y_m.$ $$H_0: \mathbb{E}X = \mathbb{E}Y,$$ т.е. матожидания выборок равны.
Альтернативные гипотезы: $$H_1: \mathbb{E}X \neq \mathbb{E}Y \quad (\text{двухсторонняя гипотеза})$$ $$H_1': \mathbb{E}X > \mathbb{E}Y \quad (\text{односторонняя гипотеза})$$ $$H_1'': \mathbb{E}X < \mathbb{E}Y \quad (\text{односторонняя гипотеза})$$
Рассмотрим, как отличаются средние рейтинги комедий и драм (считаем выборки независимыми, так как у них всего 2 пересечения).
Python:
print('mean(X) =', statistics.mean(comedy_df.average_rating),
'vs mean(Y) =', statistics.mean(drama_df.average_rating))
_, pvalue = mannwhitneyu(comedy_df.average_rating, drama_df.average_rating, alternative = 'two-sided')
print('E(X) != E(Y): pvalue =', pvalue)
_, pvalue = mannwhitneyu(comedy_df.average_rating, drama_df.average_rating, alternative = 'greater')
print('E(X) > E(Y): pvalue =', pvalue)
_, pvalue = mannwhitneyu(comedy_df.average_rating, drama_df.average_rating, alternative = 'less')
print('E(X) < E(Y): pvalue =', pvalue)
mean(X) = 4.024275862068966 vs mean(Y) = 4.029975369458128 E(X) != E(Y): pvalue = 0.38714692405450324 E(X) > E(Y): pvalue = 0.8065933038563693 E(X) < E(Y): pvalue = 0.19357346202725162
R:
%%R
cat("E(X) =", mean(comedy_df$average_rating),
"vs E(Y) =", mean(mean(drama_df$average_rating)))
cat("\nE(X) == E(Y)\n")
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "two.sided", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "two.sided", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "two.sided", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
cat("\nE(X) < E(Y)\n")
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "greater", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "greater", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "greater", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
cat("\nE(X) > E(Y)\n")
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "less", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "less", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "less", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
E(X) = 4.024276 vs E(Y) = 4.029975
E(X) == E(Y)
0.90: pvalue = 0.3871469
0.95: pvalue = 0.3871469
0.99: pvalue = 0.3871469
E(X) < E(Y)
0.90: pvalue = 0.8065933
0.95: pvalue = 0.8065933
0.99: pvalue = 0.8065933
E(X) > E(Y)
0.90: pvalue = 0.1935735
0.95: pvalue = 0.1935735
0.99: pvalue = 0.1935735
Средний рейтинг комедий меньше среднего рейтинга драм (более "серьёзные" произведения оцениваются выше).
Даны выборки $X_1, \dots, X_n \sim N(\mu_1, \sigma_1)$ и $Y_1, \dots, Y_m \sim N(\mu_2, \sigma_2).$ $$H_0: DX = DY,$$ т.е. дисперсии выборок равны.
Альтернативные гипотезы: $$H_1: DX \neq DY \quad (\text{двухсторонняя гипотеза})$$ $$H_1': DX > DY \quad (\text{односторонняя гипотеза})$$ $$H_1'': DX < DY \quad (\text{односторонняя гипотеза})$$
Рассмотрим, как отличаются разброс возрастов авторов-женщин и авторов-мужчин.
Python:
На Python нет реализации критерия Фишера, поэтому сделаем её сами.
def f_test(x, y, alternative='two_sided'):
df1 = len(x) - 1
df2 = len(y) - 1
f = x.var() / y.var()
if alternative == 'greater':
p_val = 1.0 - stats.f.cdf(f, df1, df2)
elif alternative == 'less':
p_val = stats.f.cdf(f, df1, df2)
else:
p_val = 2.0 * (1.0 - stats.f.cdf(f, df1, df2))
print(' pvalue = ', p_val)
print('var(X) =', statistics.variance(df1[df1.gender == 'M'].age),
'vs var(Y) =', statistics.variance(df1[df1.gender == 'F'].age))
print('D(X) != D(Y):', end='')
f_test(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, alternative='two-sided')
print('D(X) > D(Y):', end='')
f_test(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, alternative='greater')
print('D(X) < D(Y):', end='')
f_test(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, alternative='less')
var(X) = 119.11199528440908 vs var(Y) = 118.34482758620689 D(X) != D(Y): pvalue = 0.9983269033717543 D(X) > D(Y): pvalue = 0.49916345168587717 D(X) < D(Y): pvalue = 0.5008365483141228
R:
%%R
cat("var(X) =", var(df1[df1$gender %like% 'M', "age"]),
"vs var(Y) =", var(df1[df1$gender %like% 'F', "age"]))
cat("\nD(X) == D(Y)\n")
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "two.sided", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "two.sided", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "two.sided", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
cat("\nD(X) > D(Y)\n")
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "greater", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "greater", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "greater", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
cat("\nD(X) < D(Y)\n")
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "less", conf.level = 0.9)
cat(" 0.90: pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "less", conf.level = 0.95)
cat(" 0.95: pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "less", conf.level = 0.99)
cat(" 0.99: pvalue =", res$p.value, '\n')
var(X) = 119.112 vs var(Y) = 118.3448
D(X) == D(Y)
0.90: pvalue = 0.9983269
0.95: pvalue = 0.9983269
0.99: pvalue = 0.9983269
D(X) > D(Y)
0.90: pvalue = 0.4991635
0.95: pvalue = 0.4991635
0.99: pvalue = 0.4991635
D(X) < D(Y)
0.90: pvalue = 0.5008365
0.95: pvalue = 0.5008365
0.99: pvalue = 0.5008365
Дисперсия возрастов у мужчин и женщин равны.
Даны выборки $X_1, \dots, X_n$ и $Y_1, \dots, Y_m.$ $$H_0: DX = DY,$$ т.е. дисперсии выборок равны.
Альтернативные гипотезы: $$H_1: DX \neq DY \quad (\text{двухсторонняя гипотеза})$$
Рассмотрим, как отличаются дисперсии оценок комедий и драм.
Python:
print('var(X) =', np.var(comedy_df.average_rating),
'vs var(Y) =', np.var(drama_df.average_rating))
_, pvalue = levene(comedy_df.average_rating, drama_df.average_rating)
print('D(X) != D(Y): pvalue =', pvalue)
var(X) = 0.07031688941736028 vs var(Y) = 0.08652389101895216 D(X) != D(Y): pvalue = 0.01462297954409574
R:
%%R
xor_copy <- df2[xor(df2$genres %like% "Comedy", df2$genres %like% "Drama"), ]
xor_copy[xor_copy$genres %like% "Comedy", "is_comedy"] <- 1
xor_copy[xor_copy$genres %like% "Drama", "is_comedy"] <- 0
and_copy <- df2[rep(which(df2$genres %like% "Comedy" & df2$genres %like% "Drama"), 2), ]
and_copy[c(1, 2), "is_comedy"] <- 1
and_copy[c(3, 4), "is_comedy"] <- 0
res_copy <- rbind(xor_copy, and_copy)
%%R
cat("var(X) =", var(comedy_df$average_rating),
"vs var(Y) =", var(drama_df$average_rating), '\n\n')
leveneTest(res_copy$average_rating, res_copy$is_comedy)
var(X) = 0.0708052 vs var(Y) = 0.08673753
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 5.9993 0.01462 *
549
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Дисперсии оценок драм и комедий не равны.
Даны выборки $X_1, \dots, X_n \sim N(\mu_1, \sigma_1)$ и $Y_1, \dots, Y_m \sim N(\mu_2, \sigma_3).$ $$H_0: DX = DY,$$ т.е. дисперсии выборок равны.
Альтернативные гипотезы: $$H_1: DX \neq DY \quad (\text{двухсторонняя гипотеза})$$ $$H_1': DX > DY \quad (\text{односторонняя гипотеза})$$ $$H_1'': DX < DY \quad (\text{односторонняя гипотеза})$$
Рассмотрим, как отличаются отклонения возрастов авторов-женщин и авторов-мужчин.
Python:
print('var(X) =', np.var(df1[df1.gender == 'M'].age),
'vs var(Y) =', np.var(df1[df1.gender == 'F'].age))
_, pvalue = bartlett(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age)
print('D(X) != D(Y): pvalue =', pvalue)
var(X) = 118.09394404266199 vs var(Y) = 116.30439952437574 D(X) != D(Y): pvalue = 0.9775499143207277
R:
%%R
cat("var(X) =", var(df1[df1$gender %like% "F", "age"]),
"vs var(Y) =", var(df1[df1$gender %like% "M", "age"]))
res <- bartlett.test(df1$age ~ df1$gender)
cat("\nD(X) != D(Y): pvalue =", res$p.value)
var(X) = 118.3448 vs var(Y) = 119.112 D(X) != D(Y): pvalue = 0.9775499
Дисперсии возрастов у мужчин и женщин одинаковые.
Даны выборки $X_1, \dots, X_n$ и $Y_1, \dots, Y_m.$ $$H_0: DX = DY,$$ т.е. дисперсии выборок равны.
Альтернативные гипотезы: $$H_1: DX \neq DY \quad (\text{двухсторонняя гипотеза})$$ $$H_1': DX > DY \quad (\text{односторонняя гипотеза})$$ $$H_1'': DX < DY \quad (\text{односторонняя гипотеза})$$
Рассмотрим, как отличаются отклонения оценок комедий и драм.
Python:
print('var(X) =', np.var(comedy_df.average_rating),
'vs var(Y) =', np.var(drama_df.average_rating))
_, pvalue = fligner(comedy_df.average_rating, drama_df.average_rating)
print('D(X) != D(Y): pvalue =', pvalue)
var(X) = 0.07031688941736028 vs var(Y) = 0.08652389101895216 D(X) != D(Y): pvalue = 0.012561899939530143
R:
%%R
cat("var(X) =", var(comedy_df$average_rating),
"vs var(Y) =", var(drama_df$average_rating), '\n\n')
res <- fligner.test(res_copy$average_rating, res_copy$is_comedy)
cat("\nD(X) != D(Y): pvalue =", res$p.value)
var(X) = 0.0708052 vs var(Y) = 0.08673753 D(X) != D(Y): pvalue = 0.0125619
Дисперсии оценок драм и комедий не равны.
Рассмотрим, как коррелируют количественные данные в обоих датасетах.
Python:
print("df1:")
print("birth_year x age:")
print(" ", pearsonr(df1.birth_year, df1.age))
print("birth_year x first_published:")
print(" ", pearsonr(df1.birth_year, df1.first_published))
print("birth_year x sales:")
print(" ", pearsonr(df1.birth_year, df1.sales))
print("age x first_published:")
print(" ", pearsonr(df1.age, df1.first_published))
print("age x sales:")
print(" ", pearsonr(df1.age, df1.sales))
print("first_published x sales:")
print(" ", pearsonr(df1.first_published, df1.sales))
print("df2:")
print("average_rating x number_of_ratings:")
print(" ", pearsonr(df2.average_rating, df2.number_of_ratings))
df1:
birth_year x age:
PearsonRResult(statistic=-0.16111193049698522, pvalue=0.033173980197517605)
birth_year x first_published:
PearsonRResult(statistic=0.9858049406054875, pvalue=5.245782166441118e-136)
birth_year x sales:
PearsonRResult(statistic=-0.1273678745473681, pvalue=0.093017306705826)
age x first_published:
PearsonRResult(statistic=0.006876377241550452, pvalue=0.928036833440516)
age x sales:
PearsonRResult(statistic=0.02408953122716002, pvalue=0.7516710557994336)
first_published x sales:
PearsonRResult(statistic=-0.12495272423683916, pvalue=0.0994374716756319)
df2:
average_rating x number_of_ratings:
PearsonRResult(statistic=0.028439105170956217, pvalue=0.004453284358105489)
R:
%%R
cat("df1:\n")
cat("birth_year x age:\n ")
print(cor.test(df1$birth_year, df1$age, method="pearson"))
cat("birth_year x first_published:\n ")
print(cor.test(df1$birth_year, df1$first_published, method="pearson"))
cat("birth_year x sales:\n ")
print(cor.test(df1$birth_year, df1$sales, method="pearson"))
cat("age x first_published:\n ")
print(cor.test(df1$age, df1$first_published, method="pearson"))
cat("age x sales:\n ")
print(cor.test(df1$age, df1$sales, method="pearson"))
cat("first_published x sales:\n ")
print(cor.test(df1$first_published, df1$sales, method="pearson"))
cat("\ndf2:\n")
cat("average_rating x number_of_ratings:\n ")
print(cor.test(df2$average_rating, df2$number_of_ratings, method="pearson"))
df1:
birth_year x age:
Pearson's product-moment correlation
data: df1$birth_year and df1$age
t = -2.1471, df = 173, p-value = 0.03317
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.30223176 -0.01308145
sample estimates:
cor
-0.1611119
birth_year x first_published:
Pearson's product-moment correlation
data: df1$birth_year and df1$first_published
t = 77.228, df = 173, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.9809071 0.9894530
sample estimates:
cor
0.9858049
birth_year x sales:
Pearson's product-moment correlation
data: df1$birth_year and df1$sales
t = -1.689, df = 173, p-value = 0.09302
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.27059828 0.02137919
sample estimates:
cor
-0.1273679
age x first_published:
Pearson's product-moment correlation
data: df1$age and df1$first_published
t = 0.090447, df = 173, p-value = 0.928
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1416112 0.1550613
sample estimates:
cor
0.006876377
age x sales:
Pearson's product-moment correlation
data: df1$age and df1$sales
t = 0.31694, df = 173, p-value = 0.7517
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1246992 0.1718187
sample estimates:
cor
0.02408953
first_published x sales:
Pearson's product-moment correlation
data: df1$first_published and df1$sales
t = -1.6565, df = 173, p-value = 0.09944
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.26832226 0.02383215
sample estimates:
cor
-0.1249527
df2:
average_rating x number_of_ratings:
Pearson's product-moment correlation
data: df2$average_rating and df2$number_of_ratings
t = 2.8448, df = 9998, p-value = 0.004453
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.008843965 0.048012413
sample estimates:
cor
0.02843911
Python:
print("df1:")
print("birth_year x age:")
print(" ", spearmanr(df1.birth_year, df1.age))
print("birth_year x first_published:")
print(" ", spearmanr(df1.birth_year, df1.first_published))
print("birth_year x sales:")
print(" ", spearmanr(df1.birth_year, df1.sales))
print("age x first_published:")
print(" ", spearmanr(df1.age, df1.first_published))
print("age x sales:")
print(" ", spearmanr(df1.age, df1.sales))
print("first_published x sales:")
print(" ", spearmanr(df1.first_published, df1.sales))
print("df2:")
print("average_rating x number_of_ratings:")
print(" ", spearmanr(df2.average_rating, df2.number_of_ratings))
df1:
birth_year x age:
SignificanceResult(statistic=-0.22161108578341485, pvalue=0.0032049243336212097)
birth_year x first_published:
SignificanceResult(statistic=0.9286580663365899, pvalue=2.0088541923201725e-76)
birth_year x sales:
SignificanceResult(statistic=-0.16708646133214528, pvalue=0.02710059167888117)
age x first_published:
SignificanceResult(statistic=0.11232645210830353, pvalue=0.1388790585617484)
age x sales:
SignificanceResult(statistic=0.025341821629909604, pvalue=0.7392159748274548)
first_published x sales:
SignificanceResult(statistic=-0.1923388833606713, pvalue=0.010771516916827907)
df2:
average_rating x number_of_ratings:
SignificanceResult(statistic=-0.15452971557803954, pvalue=1.733872942153623e-54)
R:
%%R
cat("df1:\n")
cat("birth_year x age:\n ")
print(cor.test(df1$birth_year, df1$age, method="spearman"))
cat("birth_year x first_published:\n ")
print(cor.test(df1$birth_year, df1$first_published, method="spearman"))
cat("birth_year x sales:\n ")
print(cor.test(df1$birth_year, df1$sales, method="spearman"))
cat("age x first_published:\n ")
print(cor.test(df1$age, df1$first_published, method="spearman"))
cat("age x sales:\n ")
print(cor.test(df1$age, df1$sales, method="spearman"))
cat("first_published x sales:\n ")
print(cor.test(df1$first_published, df1$sales, method="spearman"))
cat("\ndf2:\n")
cat("average_rating x number_of_ratings:\n ")
print(cor.test(df2$average_rating, df2$number_of_ratings, method="spearman"))
df1:
birth_year x age:
Spearman's rank correlation rho
data: df1$birth_year and df1$age
S = 1091143, p-value = 0.003205
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.2216111
birth_year x first_published:
Spearman's rank correlation rho
data: df1$birth_year and df1$first_published
S = 63723, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.9286581
birth_year x sales:
Spearman's rank correlation rho
data: df1$birth_year and df1$sales
S = 1042442, p-value = 0.0271
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.1670865
age x first_published:
Spearman's rank correlation rho
data: df1$age and df1$first_published
S = 792870, p-value = 0.1389
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.1123265
age x sales:
Spearman's rank correlation rho
data: df1$age and df1$sales
S = 870565, p-value = 0.7392
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.02534182
first_published x sales:
Spearman's rank correlation rho
data: df1$first_published and df1$sales
S = 1064997, p-value = 0.01077
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.1923389
df2:
average_rating x number_of_ratings:
Spearman's rank correlation rho
data: df2$average_rating and df2$number_of_ratings
S = 1.9242e+11, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.1545297
Python:
print("df1:")
print("birth_year x age:")
print(" ", kendalltau(df1.birth_year, df1.age))
print("birth_year x first_published:")
print(" ", kendalltau(df1.birth_year, df1.first_published))
print("birth_year x sales:")
print(" ", kendalltau(df1.birth_year, df1.sales))
print("age x first_published:")
print(" ", kendalltau(df1.age, df1.first_published))
print("age x sales:")
print(" ", kendalltau(df1.age, df1.sales))
print("first_published x sales:")
print(" ", kendalltau(df1.first_published, df1.sales))
print("df2:")
print("average_rating x number_of_ratings:")
print(" ", kendalltau(df2.average_rating, df2.number_of_ratings))
df1:
birth_year x age:
SignificanceResult(statistic=-0.1487202050324284, pvalue=0.004121070172677824)
birth_year x first_published:
SignificanceResult(statistic=0.7829564425934482, pvalue=1.949443685647089e-52)
birth_year x sales:
SignificanceResult(statistic=-0.11886744168201617, pvalue=0.02413018270093875)
age x first_published:
SignificanceResult(statistic=0.07978447872205262, pvalue=0.12306653247813049)
age x sales:
SignificanceResult(statistic=0.011638379456385595, pvalue=0.8264558116845667)
first_published x sales:
SignificanceResult(statistic=-0.1350582730200036, pvalue=0.010250945964161366)
df2:
average_rating x number_of_ratings:
SignificanceResult(statistic=-0.1030334892459667, pvalue=2.526856184366521e-53)
R:
%%R
cat("df1:\n")
cat("birth_year x age:\n ")
print(cor.test(df1$birth_year, df1$age, method="kendall"))
cat("birth_year x first_published:\n ")
print(cor.test(df1$birth_year, df1$first_published, method="kendall"))
cat("birth_year x sales:\n ")
print(cor.test(df1$birth_year, df1$sales, method="kendall"))
cat("age x first_published:\n ")
print(cor.test(df1$age, df1$first_published, method="kendall"))
cat("age x sales:\n ")
print(cor.test(df1$age, df1$sales, method="kendall"))
cat("first_published x sales:\n ")
print(cor.test(df1$first_published, df1$sales, method="kendall"))
cat("\ndf2:\n")
cat("average_rating x number_of_ratings:\n ")
print(cor.test(df2$average_rating, df2$number_of_ratings, method="kendall"))
df1:
birth_year x age:
Kendall's rank correlation tau
data: df1$birth_year and df1$age
z = -2.8687, p-value = 0.004121
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
-0.1487202
birth_year x first_published:
Kendall's rank correlation tau
data: df1$birth_year and df1$first_published
z = 15.239, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
0.7829564
birth_year x sales:
Kendall's rank correlation tau
data: df1$birth_year and df1$sales
z = -2.2551, p-value = 0.02413
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
-0.1188674
age x first_published:
Kendall's rank correlation tau
data: df1$age and df1$first_published
z = 1.542, p-value = 0.1231
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
0.07978448
age x sales:
Kendall's rank correlation tau
data: df1$age and df1$sales
z = 0.21925, p-value = 0.8265
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
0.01163838
first_published x sales:
Kendall's rank correlation tau
data: df1$first_published and df1$sales
z = -2.5672, p-value = 0.01025
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
-0.1350583
df2:
average_rating x number_of_ratings:
Kendall's rank correlation tau
data: df2$average_rating and df2$number_of_ratings
z = -15.372, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
-0.1030335
Видна высокая корреляция между birth_year и first_published в первом датасете.
Тест $\chi^2$ используется, для проверки наличия связи между двумя категориальными признаками ($X$ и $Y$), первый из которых принимает $k$ значений, а второй - $l$. Соответствующие частоты заносят в таблицу сопряжённости:
| $x_1$ | $x_2$ | ... | $x_l$ | |
|---|---|---|---|---|
| $y_1$ | $f_{11}$ | $f_{12}$ | ... | $f_{1l}$ |
| $y_2$ | $f_{21}$ | $f_{22}$ | ... | $f_{2l}$ |
| ... | ... | |||
| $y_k$ | $f_{k1}$ | $f_{k2}$ | ... | $f_{kl}$ |
т.е. переменные $X$ и $Y$ независимы.
Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.
Проверим, зависит ли рейтинг от жанра. Для этого разделим шкалу оценивания (0 - 5) на 5 равных частей и cоставим таблицу сопряжённости.
Python:
# делим шкалу оценивания
df2['rating_cut'] = pd.cut(df2['average_rating'], bins=[-0.1, 1, 2, 3, 4, 5], labels=[1, 2, 3, 4, 5])
# отделяем жанры друг от друга
df2_explode = df2.explode('genres')
# составляем таблицу сопряжённости
cross_df = pd.crosstab(df2_explode.genres, df2_explode.rating_cut)
print('Количество строк:', cross_df.shape[0])
cross_df.head(5)
Количество строк: 617
| rating_cut | 2 | 3 | 4 | 5 |
|---|---|---|---|---|
| genres | ||||
| 12th Century | 0 | 0 | 0 | 1 |
| 15th Century | 0 | 0 | 1 | 1 |
| 16th Century | 0 | 0 | 6 | 2 |
| 17th Century | 0 | 0 | 4 | 3 |
| 18th Century | 0 | 0 | 12 | 4 |
chi2_contingency(cross_df)
Chi2ContingencyResult(statistic=5118.839776700974, pvalue=2.868874482754211e-304, dof=1848, expected_freq=array([[6.90095406e-05, 9.14376413e-04, 4.34242534e-01, 5.64774080e-01],
[1.38019081e-04, 1.82875283e-03, 8.68485068e-01, 1.12954816e+00],
[5.52076325e-04, 7.31501130e-03, 3.47394027e+00, 4.51819264e+00],
...,
[2.76038162e-04, 3.65750565e-03, 1.73697014e+00, 2.25909632e+00],
[6.90095406e-05, 9.14376413e-04, 4.34242534e-01, 5.64774080e-01],
[2.89840070e-03, 3.84038093e-02, 1.82381864e+01, 2.37205114e+01]]))
R:
%%R
# делим шкалу оценивания
df2$rating_cut <- cut(df2$average_rating, breaks=c(-0.1, 1, 2, 3, 4, 5), labels=c(1, 2, 3, 4, 5))
# отделяем жанры друг от друга
df2_explode <- df2 %>% tidyr::unnest(genres)
# составляем таблицу сопряжённости
cross_df <- table(df2_explode$genres, df2_explode$rating_cut)
cat("Количество строк:", nrow(cross_df), "\n")
head(cross_df, n=5)
Количество строк: 617
1 2 3 4 5
12thCentury 0 0 0 0 1
15thCentury 0 0 0 1 1
16thCentury 0 0 0 6 2
17thCentury 0 0 0 4 3
18thCentury 0 0 0 12 4
%%R
chisq.test(cross_df)
Pearson's Chi-squared test data: cross_df X-squared = NaN, df = 2464, p-value = NA
Рейтинг и жанр зависимы.
Точный тест Фишера используется, для проверки наличия связи между двумя категориальными признаками ($X$ и $Y$), каждый из которых принимает по $2$ значения:
| $x_1$ | $x_2$ | |
|---|---|---|
| $y_1$ | $f_{11}$ | $f_{12}$ |
| $y_2$ | $f_{21}$ | $f_{22}$ |
т.е. переменные $X$ и $Y$ независимы.
Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.
Фишера можно применять на меньших выборках данных, чем критерий $\chi^2$. Воспользуемся этим, чтобы найти пересечение используемых датасетов и проверить, существует ли связь между оценкой, которую получило произведение (больше или меньше медианы), и его продажами (больше или меньше медианы).
Python:
ultra_df = pd.merge(df1, df2, how='inner', on=['title']).drop(['author_y', 'genres_x', 'rating_cut'], axis= 1).rename(columns = {'author_x':'author', 'genres_y':'genres'})
# делим шкалы оценивания
ultra_df['rating_cut'] = pd.cut(ultra_df['average_rating'],
bins=[min(ultra_df.average_rating) - 1, statistics.median(ultra_df.average_rating),
max(ultra_df.average_rating)],
labels=['rating: less', 'rating: greater'])
ultra_df['sales_cut'] = pd.cut(ultra_df['sales'],
bins=[min(ultra_df.sales) - 1, statistics.median(ultra_df.sales), max(ultra_df.sales)],
labels=['sales: less', 'sales: greater'])
# составляем таблицу сопряжённости
cross_df = pd.crosstab(ultra_df.sales_cut, ultra_df.rating_cut)
cross_df
| rating_cut | rating: less | rating: greater |
|---|---|---|
| sales_cut | ||
| sales: less | 18 | 19 |
| sales: greater | 18 | 14 |
fisher_exact(cross_df)
SignificanceResult(statistic=0.7368421052631579, pvalue=0.6308158496561179)
R:
%%R
ultra_df <- subset(merge(df1, df2, by = "title"), select = -c(genres.x, author.y, rating_cut))
names(ultra_df)[names(ultra_df) == "author.x"] <- "author"
names(ultra_df)[names(ultra_df) == "genres.y"] <- "genres"
# делим шкалы оценивания
ultra_df$rating_cut <- cut(ultra_df$average_rating,
breaks=c(min(ultra_df$average_rating) - 1, median(ultra_df$average_rating),
max(ultra_df$average_rating)),
labels=c("rating: less", "rating: greater"))
ultra_df$sales_cut <- cut(ultra_df$sales,
breaks=c(min(ultra_df$sales) - 1, median(ultra_df$sales),
max(ultra_df$sales)),
labels=c("sales: less", "sales: greater"))
# составляем таблицу сопряжённости
cross_df <- table(ultra_df$sales_cut, ultra_df$rating_cut)
cross_df
rating: less rating: greater
sales: less 18 19
sales: greater 18 14
%%R
fisher.test(cross_df)
Fisher's Exact Test for Count Data data: cross_df p-value = 0.6308 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.2558414 2.1118582 sample estimates: odds ratio 0.7401323
Переменые независимы, связи между рейтингом и количеством проданных копий нет.
Когда независимости наблюдений нет (признак проверяется на одних и тех же субъектах), для анализа таблиц напряжённости 2x2 используется критерий МакНемара.
Рассмотрим $n$ субъектов, для каждого из которых проведено 2 теста с 2 исходами ("+" и "-"). Результаты тестов занесены в таблицу сопряжённости:
| $\text{Тест 1}$ | $\text{ВСЕГО}$ | |||
|---|---|---|---|---|
| |
|
|||
| $\text{Тест 2}$ | $+$ | $a$ | $b$ | $a + b$ |
| $-$ | $c$ | $d$ | $c + d$ | |
| $\text{ВСЕГО}$ | $a + c$ | $b + d$ | $a + b + c + d = n$ | |
т.е. маргинальные распределения для всех исходов совпадают.
Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.
Рассмотрим, изменилось ли со временем влияние пола автора на жанр произведения. Для этого найдём количество жанров, в которых большинство произведений написаны мужчиной/женщиной до/после 1990.
Python:
df1_explode = df1.explode('genres')
# до 1990
fem_genres_freq = df1_explode[(df1_explode.gender == 'F') & (df1_explode.first_published < 1990)].genres.value_counts().reset_index()
fem_genres_freq.columns = ['genre', 'fem_freq_before']
male_genres_freq = df1_explode[(df1_explode.gender == 'M') & (df1_explode.first_published < 1990)].genres.value_counts().reset_index()
male_genres_freq.columns = ['genre', 'male_freq_before']
genres_freq_before = pd.merge(fem_genres_freq, male_genres_freq, 'outer')
# после 1990
fem_genres_freq = df1_explode[(df1_explode.gender == 'F') & (df1_explode.first_published >= 1990)].genres.value_counts().reset_index()
fem_genres_freq.columns = ['genre', 'fem_freq_after']
male_genres_freq = df1_explode[(df1_explode.gender == 'M') & (df1_explode.first_published >= 1990)].genres.value_counts().reset_index()
male_genres_freq.columns = ['genre', 'male_freq_after']
genres_freq_after = pd.merge(fem_genres_freq, male_genres_freq, 'outer')
# объединяем и суммируем
genres_freq = pd.merge(genres_freq_before, genres_freq_after, 'outer').fillna(0)
genres_freq['fem_win_before'] = np.where(genres_freq.fem_freq_before >= genres_freq.male_freq_before, 1, 0)
genres_freq['fem_win_after'] = np.where(genres_freq.fem_freq_after >= genres_freq.male_freq_after, 1, 0)
genres_freq['win_win'] = np.where((genres_freq.fem_win_before == 1) & (genres_freq.fem_win_after == 1), 1, 0)
genres_freq['win_lose'] = np.where((genres_freq.fem_win_before == 1) & (genres_freq.fem_win_after == 0), 1, 0)
genres_freq['lose_win'] = np.where((genres_freq.fem_win_before == 0) & (genres_freq.fem_win_after == 1), 1, 0)
genres_freq['lose_lose'] = np.where((genres_freq.fem_win_before == 0) & (genres_freq.fem_win_after == 0), 1, 0)
genres_freq.head(n=5)
| genre | fem_freq_before | male_freq_before | fem_freq_after | male_freq_after | fem_win_before | fem_win_after | win_win | win_lose | lose_win | lose_lose | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Fiction | 23.0 | 67.0 | 23.0 | 20.0 | 0 | 1 | 0 | 0 | 1 | 0 |
| 1 | Classics | 22.0 | 68.0 | 4.0 | 10.0 | 0 | 0 | 0 | 0 | 0 | 1 |
| 2 | Childrens | 10.0 | 17.0 | 8.0 | 1.0 | 0 | 1 | 0 | 0 | 1 | 0 |
| 3 | Historical Fiction | 10.0 | 25.0 | 2.0 | 9.0 | 0 | 0 | 0 | 0 | 0 | 1 |
| 4 | Romance | 8.0 | 5.0 | 9.0 | 3.0 | 1 | 1 | 1 | 0 | 0 | 0 |
data = []
data.append([genres_freq['win_win'].sum(), genres_freq.shape[0] - genres_freq['win_lose'].sum()])
data.append([genres_freq['lose_win'].sum(), genres_freq.shape[0] - genres_freq['lose_lose'].sum()])
data
[[41, 153], [79, 138]]
print(mcnemar(data, exact=False))
pvalue 1.6456406755618912e-06 statistic 22.969827586206897
R:
%%R
df1_explode <- df1 %>% tidyr::unnest(genres)
# до 1990
fem_genres_freq <- as.data.frame(table(df1_explode[(df1_explode$gender %like% 'F') & (df1_explode$first_published < 1990),]$genres))
names(fem_genres_freq) <- c("genre", "fem_freq_before")
male_genres_freq <- as.data.frame(table(df1_explode[(df1_explode$gender %like% 'M') & (df1_explode$first_published < 1990),]$genres))
names(male_genres_freq) <- c("genre", "male_freq_before")
genres_freq_before <- merge(fem_genres_freq, male_genres_freq, all = TRUE)
# после 1990
fem_genres_freq <- as.data.frame(table(df1_explode[(df1_explode$gender %like% 'F') & (df1_explode$first_published >= 1990),]$genres))
names(fem_genres_freq) <- c("genre", "fem_freq_after")
male_genres_freq <- as.data.frame(table(df1_explode[(df1_explode$gender %like% 'M') & (df1_explode$first_published >= 1990),]$genres))
names(male_genres_freq) <- c("genre", "male_freq_after")
genres_freq_after <- merge(fem_genres_freq, male_genres_freq, all = TRUE)
# объединяем и суммируем
genres_freq <- merge(genres_freq_before, genres_freq_after, all = TRUE)
genres_freq[is.na(genres_freq)] <- 0
genres_freq$fem_win_before <- ifelse(genres_freq$fem_freq_before >= genres_freq$male_freq_before, 1, 0)
genres_freq$fem_win_after <- ifelse(genres_freq$fem_freq_after >= genres_freq$male_freq_after, 1, 0)
genres_freq$win_win <- ifelse(genres_freq$fem_win_before == 1 & genres_freq$fem_win_after == 1, 1, 0)
genres_freq$win_lose <- ifelse(genres_freq$fem_win_before == 1 & genres_freq$fem_win_after == 0, 1, 0)
genres_freq$lose_win <- ifelse(genres_freq$fem_win_before == 0 & genres_freq$fem_win_after == 1, 1, 0)
genres_freq$lose_lose <- ifelse(genres_freq$fem_win_before == 0 & genres_freq$fem_win_after == 0, 1, 0)
head(genres_freq, n=5)
genre fem_freq_before male_freq_before fem_freq_after 1 Adult 4 0 6 2 Adventure 1 15 6 3 Animals 4 4 1 4 Asia 1 2 1 5 AsianLiterature 1 1 0 male_freq_after fem_win_before fem_win_after win_win win_lose lose_win 1 4 1 1 1 0 0 2 5 0 1 0 0 1 3 4 1 0 0 1 0 4 1 0 1 0 0 1 5 0 1 1 1 0 0 lose_lose 1 0 2 0 3 0 4 0 5 0
%%R
data <- list()
data[[1]] <- c(sum(genres_freq$win_win), nrow(genres_freq) - sum(genres_freq$win_lose))
data[[2]] <- c(sum(genres_freq$lose_win), nrow(genres_freq) - sum(genres_freq$lose_lose))
%%R
mcnemar.test(matrix(unlist(data), ncol=2, byrow=TRUE))
McNemar's Chi-squared test with continuity correction data: matrix(unlist(data), ncol = 2, byrow = TRUE) McNemar's chi-squared = 22.97, df = 1, p-value = 1.646e-06
Со временем количество женров, в которых пишут преимущественно женщины, изменилось.
Если наблюдения сгруппированы по $k$ стратам, то используется тест Кохрана-Мантеля-Хензеля. В таком случае рассматривают $k$ таблиц сопряжённости 2x2:
| $x_1$ | $x_2$ | $\text{ВСЕГО}$ | |
|---|---|---|---|
| $y_1$ | $a_i$ | $b_i$ | $a_i+b_i$ |
| $y_2$ | $c_i$ | $d_i$ | $c_i+d_i$ |
| $\text{ВСЕГО}$ | $a_i+c_i$ | $b_i+d_i$ | $t_i=a_i+b_i+c_i+d_i$ |
т.е. шансы равны (между признаками нет связи).
Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.
Проверим, как связаны пол автора и жанр литературы (художестенная/нехудожественная). При этом поделим наблюдения на 3 части по возрасту автора: молодёжь (до 30 лет), взрослые (от 31 до 55 лет) и пожилые (старше 55 лет).
Python:
df1_explode['age_cut'] = pd.cut(df1_explode['age'],
bins=[min(df1_explode.age) - 1, 30,
55, max(df1_explode.age)],
labels=['young', 'adult', 'senior'])
result = CMH(df1_explode[(df1_explode.genres == 'Nonfiction') |
(df1_explode.genres == 'Fiction')],
'gender', 'genres', stratifier='age_cut')
print(result)
Cochran-Mantel-Haenszel Chi2 test
"gender" x "genres", stratified by "age_cut"
Cochran-Mantel-Haenszel M^2 = 0.2914, dof = 1, p-value = 0.5893
R:
%%R
df1_explode$age_cut <- cut(df1_explode$age,
breaks=c(min(df1_explode$age) - 1, 30,
55, max(df1_explode$age)),
labels=c("young", "adult", "senior"))
mantelhaen_array <- array(
c(nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "young", ]),
nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "young", ]),
nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "young", ]),
nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "young", ]),
nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "adult", ]),
nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "adult", ]),
nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "adult", ]),
nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "adult", ]),
nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "senior", ]),
nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "senior", ]),
nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "senior", ]),
nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "senior", ])),
dim = c(2, 2, 3),
dimnames = list(
gender = c("F", "M"),
genres = c("Fiction", "Nonfiction"),
age_cut = c("young", "adult", "senior")))
%%R
mantelhaen.test(mantelhaen_array)
Mantel-Haenszel chi-squared test with continuity correction
data: mantelhaen_array
Mantel-Haenszel X-squared = 0.11742, df = 1, p-value = 0.7318
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
0.5676904 2.6731700
sample estimates:
common odds ratio
1.231882
Между полом автора и выбором жанра нет связи.
Простейший способ графически изобразить попарную корреляцию между величинами - корреляционная матрица. Построим такие матрицы для обоих датасетов.
Python:
data1 = df1[['birth_year', 'age', 'first_published', 'sales']].dropna()
fig = plt.figure(figsize=(7,6))
ax = sns.heatmap(data1.corr(),
vmin=-1, vmax=1, cmap='RdBu',
annot = True,
annot_kws={"fontsize":9})
plt.xticks(fontsize=10)
plt.yticks(fontsize=10)
cbar = ax.collections[0].colorbar
cbar.ax.tick_params(labelsize=10)
plt.title('Корреляционная матрица: df1', fontsize=10, fontweight='bold')
plt.show()
data2 = df2[['average_rating', 'number_of_ratings']].dropna()
fig = plt.figure(figsize=(5,4))
ax = sns.heatmap(data2.corr(),
vmin=-1, vmax=1, cmap='RdBu',
annot = True,
annot_kws={"fontsize":9})
plt.xticks(fontsize=10)
plt.yticks(fontsize=10)
cbar = ax.collections[0].colorbar
cbar.ax.tick_params(labelsize=10)
plt.title('Корреляционная матрица: df2', fontsize=10, fontweight='bold')
plt.show()
R:
%%R -w 600 -h 600
data1 <- dropNA(df1[c("birth_year", "age", "first_published", "sales")])
corrplot(cor(data1),
method = "color",
tl.col = "black", addCoef.col = "gray45")
%%R -w 500 -h 500
data2 <- dropNA(df2[c('average_rating', 'number_of_ratings')])
corrplot(cor(data2),
method = "color",
tl.col = "black", addCoef.col = "gray45")
В первом датасете, ожидаемо, нашлась высокая корреляция между годом первой публикации произведения и годом рождения его автора.
Одним из способов обнаружения мультиколлинеарности является использование коэффициента инфляции дисперсии (VIF). Найдём его для обоих датасетов.
Python:
print("df1:")
X = add_constant(data1)
ds = pd.Series([variance_inflation_factor(X.values, i) for i in range(X.shape[1])],
index=X.columns)
print(ds)
df1: const 981.311602 birth_year inf age inf first_published inf sales 1.016504 dtype: float64
print("df2:")
X = add_constant(data2)
ds = pd.Series([variance_inflation_factor(X.values, i) for i in range(X.shape[1])],
index=X.columns)
print(ds)
df2: const 148.205952 average_rating 1.000809 number_of_ratings 1.000809 dtype: float64
R:
%%R
cat("df1:\n")
data1$const <- 1
# так как age и birth_year - коэффициенты-псевдонимы, используем только один из них за раз
mod <- lm(const ~ age + first_published + sales, data = data1)
cat("\n без birth_year:\n")
print(vif(mod))
mod <- lm(const ~ birth_year + first_published + sales, data = data1)
cat("\n без age:\n")
print(vif(mod))
df1:
без birth_year:
age first_published sales
1.000680 1.015962 1.016504
без age:
birth_year first_published sales
35.497757 35.475782 1.016504
%%R
cat("df2:\n")
data2$const <- 1
mod <- lm(const ~ average_rating + number_of_ratings, data = data2)
print(vif(mod))
df2:
average_rating number_of_ratings
1.000809 1.000809
VIF в первом датасете указывает на связь между birth_year, age и first_published (действительно, age = first_published - birth_year). Полученные с помощью VIF данные подтверждают то, что заметили с помощью матриц корреляции и коэффициентов корреляции.
Чтобы исследовать влияние одной или нескольких переменных категориальных переменных (факторов) на одну зависимую количественную переменную (отклик), используется дисперсионный анализ (ANOVA).
Важно, чтобы дисперсии в группах, выделенных в связи с уровнями фактора, были равны и распределение было нормальным (или достаточно близким к нему).
Рассмотрим влияние фактора с $k$ уровнями. В соответствии с уровнями фактора выделяют $k$ подвыборок размерами $n_1, \dots, n_k$ соответственно: $ X_1 = (X_{11}, \dots, X_{1n_1});\ \dots;\ X_k = (X_{k1}, \dots, X_{kn_k}).$
$$H_0: \mathbb{E}(X_1) = \dots = \mathbb{E}(X_k),$$т.е. средние выборок равны между собой (фактор не оказывает значимого влияния).
Если полученный с помощью критерия $pvalue$ больше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.
Рассмотрим, зависит ли количество оценок от рейтинга.
Python:
group1 = df2.loc[df2.rating_cut == 1].number_of_ratings
group2 = df2.loc[df2.rating_cut == 2].number_of_ratings
group3 = df2.loc[df2.rating_cut == 3].number_of_ratings
group4 = df2.loc[df2.rating_cut == 4].number_of_ratings
group5 = df2.loc[df2.rating_cut == 5].number_of_ratings
f_oneway(group1, group2, group3, group4, group5)
F_onewayResult(statistic=3.2063572589018983, pvalue=0.012197597739233967)
R:
%%R
summary(aov(number_of_ratings ~ rating_cut, df2))
Df Sum Sq Mean Sq F value Pr(>F) rating_cut 4 1.500e+12 3.75e+11 3.206 0.0122 * Residuals 9995 1.169e+15 1.17e+11 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Рейтинг книги влияет на количество оценок. Действительно, только произведения с маленьким количеством оценок получают низкий рейтинг:
print(max(df2[df2.rating_cut == 1].number_of_ratings))
print(max(df2[df2.rating_cut == 2].number_of_ratings))
print(max(df2[df2.rating_cut == 3].number_of_ratings))
print(max(df2[df2.rating_cut == 4].number_of_ratings))
print(max(df2[df2.rating_cut == 5].number_of_ratings))
0 2302 58521 6158643 9278135
%%R
cat(max(df2[df2$rating_cut == 1, "number_of_ratings"]), "\n")
cat(max(df2[df2$rating_cut == 2, "number_of_ratings"]), "\n")
cat(max(df2[df2$rating_cut == 3, "number_of_ratings"]), "\n")
cat(max(df2[df2$rating_cut == 4, "number_of_ratings"]), "\n")
cat(max(df2[df2$rating_cut == 5, "number_of_ratings"]), "\n")
0 2302 58521 6158643 9278135
Опишем случай двухфакторного анализа (большее количество факторов аналогично).
Рассмотрим влияние факторов $A$ и $B$ с $k$ и $l$ уровнями. В соответствии с уровнями факторов выделяют $k*l$ подвыборок: $ X_{1,1}, \dots, X_{k,l}.$ В таком случае формулируют 3 нулевых гипотезы: $$H_0: \mathbb{E}(X_{1i}) = \dots = \mathbb{E}(X_{ki}), \forall i$$ т.е. средние под влиянием фактора $A$ равны (фактор $A$ не оказывает значимого влияния). $$H_0: \mathbb{E}(X_{i1}) = \dots = \mathbb{E}(X_{il}), \forall i$$ т.е. средние под влиянием фактора $B$ равны (фактор $B$ не оказывает значимого влияния).
Если полученный с помощью критерия $pvalue$ больше 0.05, то соответствующая нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.
Проверим, влияют ли язык и пол и возраст автора на продажи.
Python:
df1['age_cut'] = pd.cut(df1['age'],
bins=[min(df1.age) - 1, 30,
55, max(df1.age)],
labels=['young', 'adult', 'senior'])
model = ols('sales ~ language + gender + age_cut + language:gender + language:age_cut + gender:age_cut + language:gender:age_cut', data=df1).fit()
sm.stats.anova_lm(model, typ=2)
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 15, but rank is 4
warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 1, but rank is 0
warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 2, but rank is 0
warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 15, but rank is 4
warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 30, but rank is 6
warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 2, but rank is 0
warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 30, but rank is 16
warnings.warn('covariance of constraints does not have full '
| sum_sq | df | F | PR(>F) | |
|---|---|---|---|---|
| language | 318.481500 | 15.0 | 0.025895 | 0.998688 |
| gender | NaN | 1.0 | NaN | NaN |
| age_cut | NaN | 2.0 | NaN | NaN |
| language:gender | 7208.630522 | 15.0 | 0.586121 | 0.673192 |
| language:age_cut | 9493.975688 | 30.0 | 0.385969 | 0.887081 |
| gender:age_cut | NaN | 2.0 | NaN | NaN |
| language:gender:age_cut | 17166.811598 | 30.0 | 0.697901 | 0.792587 |
| Residual | 118069.326217 | 144.0 | NaN | NaN |
R:
%%R
df1$age_cut <- cut(df1$age,
breaks=c(min(df1$age) - 1, 30,
55, max(df1$age)),
labels=c("young", "adult", "senior"))
%%R
summary(aov(sales ~ language * gender * age_cut, df1))
Df Sum Sq Mean Sq F value Pr(>F) language 15 7179 478.6 0.584 0.884 gender 1 27 27.5 0.033 0.855 age_cut 2 2181 1090.4 1.330 0.268 language:gender 5 2733 546.6 0.667 0.649 language:age_cut 5 4838 967.6 1.180 0.322 gender:age_cut 2 299 149.6 0.183 0.833 Residuals 144 118069 819.9
Вышеперечисленные факторы и их сочетания не влияют на продажи.
Попробуем понять, как количество оценок влияет на рейтинг. Обозначим эти переменные как $y$ и $x$ соответственно. Сравним результаты регрессии методом наименьших квадратов следующих линейных по параметру моделей: линейной ($y=wx+b$), гиберболы ($y=w\frac{1}{x+1}+b$), корня ($y=w\sqrt{x}+b$) и логарифмической ($y=w\ln (x+1)+b$).
Python:
x = df2.number_of_ratings
y = df2.average_rating
x_no0 = df2[(df2.number_of_ratings != 0) & (df2.average_rating != 0)].number_of_ratings
y_no0 = df2[(df2.number_of_ratings != 0) & (df2.average_rating != 0)].average_rating
lin_model = sm.OLS(y, sm.add_constant(x)).fit()
print(lin_model.summary())
#X = np.column_stack((np.ones(x.shape[0]), x, x**2, x**3))
hyper_model = sm.OLS(y, sm.add_constant(1 / (x + 1))).fit()
print(hyper_model.summary())
sqrt_model = sm.OLS(y, sm.add_constant(x ** 0.5)).fit()
print(sqrt_model.summary())
log_model = sm.OLS(y, sm.add_constant(np.log(x + 1))).fit()
print(log_model.summary())
OLS Regression Results
==============================================================================
Dep. Variable: average_rating R-squared: 0.001
Model: OLS Adj. R-squared: 0.001
Method: Least Squares F-statistic: 8.093
Date: Mon, 11 Dec 2023 Prob (F-statistic): 0.00445
Time: 12:34:39 Log-Likelihood: -3259.3
No. Observations: 10000 AIC: 6523.
Df Residuals: 9998 BIC: 6537.
Df Model: 1
Covariance Type: nonrobust
=====================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------
const 4.0660 0.003 1170.342 0.000 4.059 4.073
number_of_ratings 2.787e-08 9.8e-09 2.845 0.004 8.67e-09 4.71e-08
==============================================================================
Omnibus: 4767.098 Durbin-Watson: 1.932
Prob(Omnibus): 0.000 Jarque-Bera (JB): 171661.269
Skew: -1.633 Prob(JB): 0.00
Kurtosis: 23.033 Cond. No. 3.67e+05
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.67e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
OLS Regression Results
==============================================================================
Dep. Variable: average_rating R-squared: 0.002
Model: OLS Adj. R-squared: 0.002
Method: Least Squares F-statistic: 19.31
Date: Mon, 11 Dec 2023 Prob (F-statistic): 1.12e-05
Time: 12:34:39 Log-Likelihood: -3253.7
No. Observations: 10000 AIC: 6511.
Df Residuals: 9998 BIC: 6526.
Df Model: 1
Covariance Type: nonrobust
=====================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------
const 4.0648 0.003 1175.209 0.000 4.058 4.072
number_of_ratings 0.2544 0.058 4.395 0.000 0.141 0.368
==============================================================================
Omnibus: 5701.374 Durbin-Watson: 1.932
Prob(Omnibus): 0.000 Jarque-Bera (JB): 272907.764
Skew: -2.042 Prob(JB): 0.00
Kurtosis: 28.264 Cond. No. 17.3
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
OLS Regression Results
==============================================================================
Dep. Variable: average_rating R-squared: 0.001
Model: OLS Adj. R-squared: 0.001
Method: Least Squares F-statistic: 6.808
Date: Mon, 11 Dec 2023 Prob (F-statistic): 0.00909
Time: 12:34:39 Log-Likelihood: -3259.9
No. Observations: 10000 AIC: 6524.
Df Residuals: 9998 BIC: 6538.
Df Model: 1
Covariance Type: nonrobust
=====================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------
const 4.0755 0.004 955.063 0.000 4.067 4.084
number_of_ratings -3.65e-05 1.4e-05 -2.609 0.009 -6.39e-05 -9.08e-06
==============================================================================
Omnibus: 4846.069 Durbin-Watson: 1.930
Prob(Omnibus): 0.000 Jarque-Bera (JB): 174716.871
Skew: -1.674 Prob(JB): 0.00
Kurtosis: 23.202 Cond. No. 388.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
OLS Regression Results
==============================================================================
Dep. Variable: average_rating R-squared: 0.030
Model: OLS Adj. R-squared: 0.030
Method: Least Squares F-statistic: 310.9
Date: Mon, 11 Dec 2023 Prob (F-statistic): 1.54e-68
Time: 12:34:39 Log-Likelihood: -3110.2
No. Observations: 10000 AIC: 6224.
Df Residuals: 9998 BIC: 6239.
Df Model: 1
Covariance Type: nonrobust
=====================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------
const 4.2239 0.009 448.912 0.000 4.205 4.242
number_of_ratings -0.0179 0.001 -17.631 0.000 -0.020 -0.016
==============================================================================
Omnibus: 5957.727 Durbin-Watson: 1.919
Prob(Omnibus): 0.000 Jarque-Bera (JB): 261509.787
Skew: -2.213 Prob(JB): 0.00
Kurtosis: 27.658 Cond. No. 26.7
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
fig, axes = plt.subplots(figsize=(10,10))
plt.scatter(x, y, alpha = 0.1)
polyline = np.linspace(min(x), max(x), max(x)-min(x))
polyline2 = np.linspace(0, 5, 1000)
axes.plot(polyline, lin_model.params[0] + lin_model.params[1] * polyline,
color="red", label="Линейная")
axes.plot(polyline, hyper_model.params[0] + hyper_model.params[1] * (1 / (polyline + 1)),
color="yellow", label="Гипербола")
axes.plot(polyline, sqrt_model.params[0] + sqrt_model.params[1] * polyline ** 0.5,
color="green", label="Корень")
axes.plot(polyline, log_model.params[0] + log_model.params[1] * np.log(polyline + 1),
color="blue", label="Логарифм")
plt.xticks(fontsize=9)
plt.yticks(fontsize=9)
#plt.xlim(-10, 0.2 * 1e6)
axes.legend(loc="best", fontsize=10)
plt.show()
R:
%%R
x <- df2$number_of_ratings
y <- df2$average_rating
lin_model <- lm(y ~ x)
print(summary(lin_model))
hyper_model <- lm(y ~ 1/(x+1))
print(summary(hyper_model))
sqrt_model <- lm(y ~ sqrt(x))
print(summary(sqrt_model))
log_model <- lm(y ~ log(x+1))
print(summary(log_model))
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-4.0660 -0.1884 0.0040 0.1928 0.9340
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.066e+00 3.474e-03 1170.342 < 2e-16 ***
x 2.787e-08 9.799e-09 2.845 0.00445 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3352 on 9998 degrees of freedom
Multiple R-squared: 0.0008088, Adjusted R-squared: 0.0007088
F-statistic: 8.093 on 1 and 9998 DF, p-value: 0.004453
Call:
lm(formula = y ~ 1/(x + 1))
Residuals:
Min 1Q Median 3Q Max
-4.0686 -0.1886 0.0114 0.1914 0.9314
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.068577 0.003354 1213 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3354 on 9999 degrees of freedom
Call:
lm(formula = y ~ sqrt(x))
Residuals:
Min 1Q Median 3Q Max
-4.0755 -0.1893 0.0075 0.1911 0.9249
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.075e+00 4.267e-03 955.063 < 2e-16 ***
sqrt(x) -3.650e-05 1.399e-05 -2.609 0.00909 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3353 on 9998 degrees of freedom
Multiple R-squared: 0.0006804, Adjusted R-squared: 0.0005805
F-statistic: 6.808 on 1 and 9998 DF, p-value: 0.009091
Call:
lm(formula = y ~ log(x + 1))
Residuals:
Min 1Q Median 3Q Max
-4.2239 -0.1767 0.0214 0.1944 0.8571
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.223915 0.009409 448.91 <2e-16 ***
log(x + 1) -0.017876 0.001014 -17.63 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3303 on 9998 degrees of freedom
Multiple R-squared: 0.03015, Adjusted R-squared: 0.03006
F-statistic: 310.9 on 1 and 9998 DF, p-value: < 2.2e-16
%%R -w 1000 -h 1000
plot(x, y)
abline(lin_model, col='red')
xs = seq(from = min(x), to = max(x), length.out = max(x)-min(x))
hyper_ys = predict(hyper_model, newdata = list(x=seq(from = min(x), to = max(x), length.out = max(x)-min(x))))
matlines(xs, hyper_ys, lwd = 2, col = 'yellow')
sqrt_ys = predict(sqrt_model, newdata = list(x=seq(from = min(x), to = max(x), length.out = max(x)-min(x))))
matlines(xs, sqrt_ys, lwd = 2, col = 'green')
log_ys = predict(log_model, newdata = list(x=seq(from = min(x), to = max(x), length.out = max(x)-min(x))))
matlines(xs, log_ys, lwd = 2, col = 'blue')
Полученные модели, особенно последние три, дают очень близкий результат (в частности, их остаточные стандартные ошибки крайне близки).
Также по коэффициентам регрессии видно, что количество оценок практически не влияет на рейтинг.